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FOREWORD 


A  library  of  subroutines  was  prepared  in  machine  language  for  the  Naval  Ordnance 
Research  Calculator.  Many  of  the  subroutines  have  been  rewritten  in  FORTRAN.  Those 
subroutines  in  the  library  which  perform  operations  on  polynomials  and  on  matrices 
have  been  documented  in  previous  reports.  Those  subroutines  which  compute  special 
functions  are  documented  herewith.  The  manuscript  for  this  report  was  completed  by 
7  November  1978. 


Released  by: 

Ralph  A.  Niemann 

Head,  Strategic  Systems  Department 


ii 


ABSTRACT 

Documentation  is  given  for  some  subroutines  which  compute  potentials  and  other 
functions.  A  set  of  subroutines  uses  rational  approximations  to  compute  Bessel  functions 
of  integral  order.  One  subroutine  uses  the  Debye  approximation  for  the  efficient 
computation  of  Bessel  functions  of  complex  argument  and  complex  order.  Empirical 
formulae  have  been  developed  to  express  the  limiting  boundaries  of  the  modes  of 
computation. 


INTRODUCTION 


On  the  Naval  Ordnance  Research  Calculator,  programs  were  coded  directly  in  machine 
language.  It  was  necessary  to  provide  subroutines  for  such  elementary  functions  as 
square  root,  sine,  cosine,  exponential,  logarithm,  and  arctangent.  With  the  advent  of 
the  FORTRAN  compiler,  versions  of  such  functions  were  available  from  the  compiler, 
and  it  was  no  longer  considered  necessary  to  provide  function  routines  with  each 
program. 

The  function  name  SORT  has  been  preempted  in  FORTRAN  for  the  square  root.  The 
function  name  CBRT  is  used  on  the  Univac  1108  computer  for  the  cube  root.  A  new 
function  routine  for  the  cube  root  has  been  prepared  for  the  CDC  6600  computer.  It 
obviates  the  inefficiency  of  exponentiation. 

Too  many  problems  have  arisen  with  the  four- quadrant  arctangent  function  routines 
of  the  Control  Data  Corporation.  Among  the  versions  and  revisions  of  the 
double-precision  DATAN2  function  routines  on  Scope  3.3.  some  have  returned  nonzero 
when  they  should  have  returned  zero,  and  some  have  returned  zero  when  they  should 
have  returned  +n.  Even  the  single-precision  ATAN2  function  routine  on  Scope  3.4 
returns  completely  erroneous  values  for  angles  outside  of  that  quadrant  which  straddles 
zero.  The  double-precision  DAT  AN  2  function  routine  on  Scope  3.4  returns  zero  when 
it  should  return  +n. 

Patches  in  FORTRAN  have  been  designed  by  A.  H.  Morris,  Jr.  to  cure  the  arctangent 
function  routines  in  Scope  3.4.  These  patches  convert  arguments  into  absolute  values 
before  series  expansion  and  restore  signs  to  the  function  after  series  expansion.  The 
patches  do  not  stand  alone,  because  they  depend  upon  the  system  function  routines 

ABS  and  AT AN. 

New  FORTRAN  function  routines  have  been  prepared  for  both  A‘rAN2  and  DATAN2. 
When  they  are  included  with  a  program  deck,  they  override  the  function  routines  with 
the  same  names  already  in  the  system.  The  new  FORTRAN  function  routines  use  the 
signs  of  arguments  directly  to  determine  the  sign  of  the  function.  They  stand  alone 
because  they  make  no  reference  to  systems  function  routines. 

The  arctangent  is  a  multiple-valued  function.  If  only  the  projections  x,y  of  a  line 
on  the  coordinate  axes  are  given,  then  the  angle  which  the  line  makes  with  the  x-axis 
is  indeterminate  modulo  2 n.  The  accepted  convention  is  to  assume  that  the  arctangent 
satisfies  the  inequality 


-  n  <  tan'*^\  §  +  7t  (1) 

Any  arctangent  function  routine  must  give  +tt  when  x  is  negative  and  y  is  zero. 
Otherwise  incorrect  results  will  be  obtained  when  real  numbers  are  included  in  a  set 
of  complex  numbers,  as  often  happens  in  mathematical  exercises. 

In  arctangent  function  routines  the  addition  theorem  is  used  for  various  centers 
of  expansion  and  the  Maclaurin  expansion  is  used  at  each  center  of  expansion.  An 
increase  in  the  number  of  centers  permits  a  decrease  in  the  number  of  terms.  In  the 
CDC  version  the  range  of  tangents  from  0  to  1  is  divided  into  sixteen  sectors.  A 
five-term  expansion  is  used  from  one  center  to  the  next  higher  center.  In  the  new 
version  the  range  of  angle  from  -|rr  to  +|n  is  divided  into  sectors  at  seven  angles 
for  which  the  tangents  are  known  to  especially  many  digits  of  accuracy.  The  series 
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expansion  runs  through  eight  terms  or  less  from  the  midpoint  between  one  pair  of 
centers  to  the  midpoint  between  the  next  higher  pair  of  centers. 

Bessel  functions  are  given  analytically  by  absolutely  convergent  series  everywhere 
in  the  complex  plane,  but  evaluation  of  the  convergent  series  by  computer  is  feasible 
only  in  a  limited  range  of  order  and  argument.  There  are  asymptotic  series  which  are 
valid  for  large  argument,  but  evaluation  of  the  asymptotic  series  also  is  feasible  only 
in  a  limited  range  of  order  and  argument.  The  Debye  asymptotic  approximation  can 
be  used  to  reduce  the  gap  in  the  range  of  order  and  argument  to  a  narrow  zone  which 
straddles  the  line  of  equal  order  and  argument.  The  zone  which  still  is  not  covered 
by  the  series  can  be  crossed  with  the  aid  of  recurrence  equations.  Evaluation  of  the 
Debye  series  requires  a  double  summation,  but  the  time  for  evaluation  is  less  than 
the  time  which  would  be  required  for  recurrence  from  the  classical  series.  Formulae 
for  the  Debye  approximation  have  been  given  by  Watson1,2  and  by  Abramowitz  and 
Stegun3,  while  explicit  recurrence  equations  for  the  terms  of  the  series  have  been 
given  by  Amos6-10. 


CUBE  ROOT 

Analysis 

One  third  of  the  exponent  of  x  is  biased  and  is  attached  to  1  to  form  an  initial 
approximation  which  is  larger  than  the  cube  root  of  x.  The  initial  approximation  is 
diminished  by  Newton-Raphson  iteration  until  the  increment  in  root  is  zero  or  positive 
from  rounding  error. 

Programming 

FUNCTION  CBRT  (X) 

********************************************************************************************* 
FORTRAN  FUNCTION  ROUTINE  FOR  CUBE  ROOT 
******************.******* *************************************************  *******  *  *********** 

The  cube  root  of  x  with  sign  is  computed  by  Newton-Raphson  iteration  and  is  stored 
in  address  CBRT. 


FOUR-QUADRANT  ARCTANGENT 

Analysis 

Let  x,  y  be  the  arguments  of  the  arctangent  function.  The  given  arguments  z,  y  can 
be  replaced  by  new  arguments  u,  v  through  the  application  of  symmetry  relations.  Let 
c  be  a  constant  which  is  added  to  the  arctangent  of  it.  v  and  let  h  be  the  center  of 
expansion  of  the  arctangent  of  u,  v. 

If  z,  y  satisfy  the  inequalities 

ll/l  §|x|  0  <  x  (2) 

then  u,  v,  c  are  given  by  the  substitutions 

u  -  +  x  V  -  +  y  c-0  (3) 
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If  x.y  satisfy  the  inequality 


(4) 


1*1  <  \y\ 


then  u,  v,  c  are  given  by  the  substitutions 

u-»+y  v  -  x  c  -»  ±  (5) 

where  the  sign  of  c  is  the  same  as  the  sign  of  y.  If  *,  y  satisfy  the  inequality 

Ivl  ^  M  X  <  0  (6) 

then  u,  v,  c  are  given  by  the  substitutions 

u  -*  x  v  ■*  y  c  ->  ±  tt  (7) 

where  the  sign  of  c  is  the  same  as  the  sign  of  y. 

The  parameters  c  and  h  are  adjusted  in  accordance  with  the  substitutions 

c  -» c  +  8  h  -» tan  8  (8) 

where  8  is  that  angle  among  the  angles 

0  ±  &  ±  g»r  ±  irr  (9) 

which  is  nearest  to  tan“l(ii/ti),  and  tan  8  is  the  corresponding  tangent  among  the 
tangents 


0 


±  (2  -  \/3) 


(10) 


Then  the  Maclaurin  expansion  is  given  by  the  equation 

,  ;  (-I)-,*-*' 

’-JLisrrrr 

for  which  the  parameter  q  is  given  by  the  equation 

v  -  hu 
9  = 


u  +  hv 


(11) 


(12) 


and  the  arctangent  function  is  given  by  the  equation 

tan-1^  j  =  c  +  tan_lq  (13) 

Multiple-precision  values  for  the  constants  c  and  h  are  to  be  found  in  the  following 
table,  which  is  derived  from  references  17  and  18. 

2  -  V3  =  0.26794  9 1 924  3 1 1 22  70647  25536  58493  87236 

=  0.57735  02691  89625  76450  91487  80502  04254 

=  0.26179  93877  99149  43653  85536  15273  29191 

|tt  =  0.52359  87755  98298  87307  71072  30546  58382 

4tt  =  0.78539  81633  97448  30961  56608  45819  87572 

|tt  =  1.57079  63267  94896  61923  13216  91639  75144 

n  =  3. 14159  26535  89793  23846  26433  83279  50288 


$ 
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Multiplication  or  division  of  the  numbers  by  small  integers  can  be  verified  by  hand 
computation. 

Programming 

FUNCTION  ATAN2  (Y,  X) 

********************************************************************************************* 
FORTRAN  FUNCTION  ROUTINE  FOR  SINGLE  PRECISION  ARCTANGENT 
********************************************************************************************* 

The  variables  x.y  are  given  in  the  arguments  X,  Y.  The  four-quadrant  arctangent 
of  y/x  is  returned  as  the  function  ATAN2. 

DOUBLE  FUNCTION  DATAN2  (Y.  X) 

********************************************* ************************************************ 
FORTRAN  FUNCTION  ROUTINE  FOR  DOUBLE  PRECISION  ARCTANGENT 
********************************************************************************************* 

The  variables  x,  y  are  given  in  the  arguments  X.Y.  The  four-quadrant  arctangent 
of  y/x  is  returned  as  the  function  DATAN2. 


POTENTIAL  OF  PLATE 


Let  a  plate  have  unit  mass  per  unit  area.  Let  x,  y,  z  be  Cartesian  coordinates  with 
origin  at  the  center  of  the  plate,  and  with  z  in  the  direction  perpendicular  to  the 
plate.  Let  r  be  the  distance  to  a  point  in  the  field  from  an  element  of  surface  ds  on 
the  plate.  The  potential  <p  at  the  field  point  is  given  by  the  equation 


^  _  J  |cfs| 

The  gradient  of  the  potential  is  given  by  the  equation 

f  Vr  |ds  | 

-v,  =  j  — 


where  Vr  is  a  unit  vector  in  the  direction  toward  the  field  point.  The  derivative  of  <p 
with  respect  to  z  is  given  by  the  equation 


dip 

dz 


Vr ds 


=  u 


(16) 


where  k  is  a  unit  vector  in  the  direction  of  increasing  z,  and  u  is  the  solid  angle  of 
the  plate.  The  derivative  of  tp  with  respect  to  z  also  is  the  potential  of  a  uniformly 
polarized  plate.  Since  the  derivative  of  a  solution  of  Laplace’s  equation  with  respect 
to  a  Cartesian  coordinate  is  itself  a  solution  of  Laplace’s  equation,  both  the  potential 
and  the  solid  angle  are  solutions  of  Laplace’s  equation. 

In  the  definition  of  a  function  at  a  point  in  the  field  it  is  convenient  to  regard  r 
as  the  distance  from  a  point  on  the  plate  to  a  point  in  the  field,  while  in  the 
transformation  of  integrals  it  is  convenient  to  regard  r  as  the  distance  from  the  point 
in  the  field  to  a  point  on  the  plate.  In  either  case  the  gradients  of  r  differ  only  in 
sign. 
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The  field  of  a  unit  current  along  the  perimeter  of  the  plate  is  defined  by  the  equation 

c£r  xVr 


-  V<p 


% 

=  O 
J 


while  for  transformation  the  field  is  expressed  by  the  equation 

-  V<fi  =  (j)  cfrxV^i  j 

Application  of  the  scalar-vector  triple  product  theorem  gives  the  equation 


fdrxV(r)  =  fdrXV(^)'I  =  ^drV(r)Xl 

sor.  Application  of  the  Stokes  1 

fdr^)xi=ldsHv(;H 


(17) 


(18) 


(19) 


where  1  is  the  identity  tensor.  Application  of  the  Stokes  theorem  gives  the  equation 

(20) 


(21) 


Application  of  the  vector-vector  triple  product  theorem  gives  the  equation 


Since  the  gradient  of  the  gradient  of  a  scalar  function  is  symmetric,  the  field  is  given 
by  the  equation 


—  V<f 


-H;) 


ds 


(22) 


Thus  the  potential  of  a  circuit  of  unit  current  is  just  the  solid  angle  of  the  circuit. 


CIRCULAR  DISK 


Analysis 

Let  a  be  the  radius  of  a  disk  of  unit  mass  per  unit  area.  Let  x,  y,  z  be  Cartesian 
coordinates  with  z  in  the  direction  of  the  axis  of  the  disk. 

The  potential  and  the  solid  angle  of  the  disk  may  be  expanded  in  a  series  of  spherical 
harmonics.  Symmetry  about  the  axis  of  the  disk  requires  that  only  symmetric  harmonics 
may  appear  in  the  series  expansion.  The  coefficients  of  the  spherical  harmonics  may 
be  derived  by  reference  to  special  series  expansions  on  the  axis  of  the  disk. 

Let  u  be  the  radial  distance  from  the  center  of  the  disk.  Then  the  potential  along 
the  axis  is  given  by  the  equation 

,-arf  -  /v=  ,  =  2nf\/az  +  z2  -  z)  (23) 

Jo  vu2  +  z2 


and  the  solid  angle  along  the  axis  is  given  by  the  equation 


d<fi 

dz 


=  27t| 


1  - 


\a2  +  z2 


(24) 


Let  r  be  the  distance  from  the  axis  of  the  disk  to  the  point  in  the  field.  Expansion  in 
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series  of  ascending  powers  of  z  and  identification  of  the  powers  of  z  with  spherical 
harmonics  lead  to  the  equations 


<?  2n\z\  2na  JL  (2n  _  1)2a«(n!)4  a  ) 

"  n?o  S^n!)2  (  a  )  P2n+1lv7^ 


(Vr2  +  2=  <  a)  (25) 


(Vr2  +  z2  <  a)  (26) 


Expansion  in  series  of  descending  powers  of  2  and  identification  of  the  powers  of  z 
with  spherical  harmonics  lead  to  the  equations 


V  M)n(2n)!  (  a  \2"+l 

V  Trant0  22nn!(n+  l)!lvr2+22/  * 

“  (~  l)n(2n+ 1 )!  j  a  \2n+2 
71  n-o  22nn!(n  +  1 )!  (  Vr^+T2  /  2" 


(vr2  + 


a)  (27) 


a)  (28) 


The  convergence  of  the  series  deteriorates  as  the  point  in  the  field  approaches  a 
sphere  of  radius  a. 

Let  a  line  be  constructed  through  the  field  point  and  perpendicular  to  the  plane 
of  the  disk.  Let  the  intersection  between  the  perpendicular  line  and  the  plane  of  the 
disk  be  the  center  of  a  circular  arc.  Let  the  reference  line  for  azimuth  angle  be  the 
extension  of  the  line  from  the  center  of  the  arc  to  the  center  of  the  disk  The  circular 
arc  intersects  the  edge  of  the  disk  at  a  point  whose  radius  from  the  center  of  the 
disk  makes  an  angle  0  with  the  reference  line.  Then  the  radius  of  the  circular  arc  is 
given  by  the  expression 

Va2  +  2ar  cos  0  +  r2  (29) 

The  derivative  of  the  radius  with  respect  to  0  is  given  by  the  expression 

-  -arSin*  -  -  (30) 

v  a2  +  2ar  cos  0  +  r2 

The  distance  of  the  circular  arc  from  the  field  point  is  given  by  the  expression 

V a2  +  2ar  cos  0  +  r2  +  22  (31) 

The  potential  of  the  disk  is  given  therefore  by  the  equation 

„  =  2ar  I"  tan  si"  * ,  ■  (32) 

Jo  \r  +  a  cos  0  /  Va2  +  2ar  cos  0  +  r2  +  c2 

The  arctangent  can  be  removed  from  the  integrand  through  an  integration  by  parts 
The  interpretation  of  the  arctangent  at  the  limits  of  integration  depends  upon  whether 
the  center  of  the  arc  is  inside  the  perimeter  of  the  disk  or  outside  the  perimeter  of 
the  disk. 

If  the  center  of  the  arc  is  inside  the  perimeter  of  the  disk,  then  the  potential  is 
given  by  the  equation 


2-nJzi  4  2a 


J.U  2c 


r  cos  ip  \ 
r  cos  0  +  r2  / 


2ar  cos  0  +  r2  +  z2  d0  (r  v  a)  (33) 
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and  the  solid  angle  is  given  by  the  equation 

+  r cos 0 


dtp 

-  —  =  ±  2rr  -  2 az 
dz 


rn  /  a  +  r  cos  0  \ 

J0  \a2  +  2ar  cos  0  +  r2/ 


dtp 


2ar  cos  tp  +  r  /  \/a2  +  2ar  cos  0  +  r2  +  z2 


(r  <  a)  (34) 


If  the  center  of  the  arc  is  outside  the  perimeter  of  the  disk,  then  the  potential  is 
given  by  the  equation 

tp  =  +  2a  f  i—s — — - — - - a2  +  2ar  cos  0  +  r2  +  z2  dtp  (a  <  r)  (35) 

J0  \aa  +  2arcos0  +  r2/  y  \  / 

and  the  solid  angle  is  given  by  the  equation 

,  +  r  cos  0 


dtp 

-  —  =  -  2az 
3z 


d0 


2ar  cos  tp  +  r2  /  Va2  +  2ar  cos  tp  +  r2  +  z^ 
A  rearrangement  of  terms  leads  to  the  equation 

dtp 


(a  <  r)  (36) 


+  2a r  cos  tp  +  r2  +  z2 


+  (a2  -  r2)z2 


Tv* 

Jo 

f 


+  2ar  cos  tp  +  r2  +  z2  dtp 


dtp 


and  to  the  equation 

dtp 
dz 


o  (a2  +  2ar  cos  0  +  r2)x/a^+~ 2ar  cos  tp  +  r2  +  z2 

r  ** 

~  2  Jo 


(a  <  r)  (37) 


’  +  2ar  cos  <f>  +  r2  +  z2 


-  (a2  -  r2)z 
The  substitution 


f 


a!0 


(a2  +  2ar  cos  tp  +  r 2)  va2  +  2ar  cos  0  +  r2  +  z2 


(a  <  r)  (38) 

(39) 


cos  0=1-2  sin2|0 

and  replacement  of  ^0  by  B  leads  to  the  expression  of  the  integrals  in  terms  of 
Legendre  elliptic  integrals.  The  first  and  second  kinds  of  elliptic  integral  are  defined 
by  the  equations 


F(tp.  k)  =  r 
Jo 


dB 


E(<p,k) 


-£ 


V 1  -  fc2sin20 

and  the  third  kind  of  elliptic  integral  is  defined  by  the  equation 

dO 


v  1  -  fc2sinzfi  dB  (40) 


11(0,  a,  fc)  =  f- - - 

Jo  ( 1  —  cx 


la  (1  -  a2sin20)v/ 1  -  fc2sin20 
If  the  modulus  a  is  defined  by  the  equation 

,  4ar 


a 


(a  +  r)2 


(41) 


(42) 
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and  the  modulus  k  is  defined  by  the  equation 


( a  +  r)*+z 2 

then  the  potential  is  given  by  the  equation 

/  2  2\ 

9=~  2  n|a|  +  2  ’■  =F(  f.  k)  +  2^(a  +  r)2  +  22£'(f,  k) 

v  (a  +  r)2  +  22 

(a  —  t)  2  2 

+  2-7— - n(f,  a.  fc) 

(a  +  r)  V(a  +  r)2  +  22 


(r  <  a)  (44) 


and  the  solid  angle  is  given  by  the  equation 


9<p  „ 

—  —  ±  2tt  — 
92 


/(a  +  r)2  +  2 


j/Xf  *)-8 


(a  -  r) 


2  +  2  2  (a  +  r)  v/(a  +  r)2  +  22  2 


II(f ,  a,  k)  (r  <  a )  (45) 


when  the  field  point  is  inside  the  perimeter,  while  the  potential  is  given  by  the  equation 

9  =  +  2-7-— ~F(f.  *)  +  2V(a  +  r)2  +  z*E(j,  k ) 
v(a  +  r)z  -f  2 2 

(a  -  r)  2 2 

,  r-v  v  /  n  /  Tf  I  \  id 


+  2^ - £  ,  =n(Sta,fe) 

(a  +  r)  V(a  +  r)2  +  22 


(a  <  r)  (46) 


and  the  solid  angle  is  given  by  the  equation 


'{ a  +  r)2  +  2 


\F{\,  k)  -  2 


(a  -  r)  2 

(a  +  r)  \/(a  +  r)2  +  22 


n(f,a,fc)  (a  <  r)  (47) 


when  the  field  point  is  outside  the  perimeter. 

The  complete  elliptic  integral  of  the  third  kind  is  expressed  in  terms  of  the  incomplete 
integrals  of  the  first  two  kinds  by  formulae  on  page  228  of  Handbook  of  Elliptic 
Integrals  by  Byrd  and  Friedman4.  If  the  angle  9  is  defined  by  the  equation 


9  =  sin’ 


a2  -  -fc2 


and  the  parameter  A(0,  Jfc)  is  defined  by  the  equation 


A (9,  -fc)  =  -  E(% ,  k)F(9,  k‘)  +  F( f ,  k)E(9,  -fc')  -  F( £,  k)F{9,  k') 

7T  |_ 


where  the  comodulus  k'  is  defined  by  the  equation 


then  the  complete  elliptic  integral  of  the  third  kind  is  given  by  the  equation 

"  a\(0,k) 


n(|.  a,  k)  = 


2  V(a2  fc2)(l  -  a2) 


When  the  moduli  are  expressed  in  terms  of  the  coordinates,  then  9  is  given  by  the 


equation 


2 


6  =  tan’1 


|a-r| 

and  the  moduli  can  be  combined  as  expressed  by  the  equation 

a  (a  +  r)  V(a  +  r)2  +  z2 


V(a* -**)(!  -  a*) 


a  -  r 


(52) 


(53) 


Thus  the  potential  and  the  solid  angle  are  reduced  to  simple  expressions  in  terms  of 
A(B.k). 

Accuracy  and  efficiency  were  determined  by  comparisons  between  computations  by 
two  formulations  on  a  common  boundary  between  their  zones  of  application. 
Optimization  of  accuracy  and  efficiency  limits  the  use  of  the  elliptic  integrals  to  an 
annular  zone  between  a  sphere  of  radius  and  a  sphere  of  radius  2a. 

Programming 

SUBROUTINE  CDSKP  (AA.  AR,  AZ.  FP) 

*************************%***********«******************************************************« 
FORTRAN  SUBROUTINE  FOR  POTENTIAL  OF  CIRCULAR  DISK 

#**^*******4^*****4r**#*#*#**************#*#*#4^********************#*******!t**«***********t* 


The  radius  a  of  the  disk  is  given  in  argument  AA,  the  distance  r  from  the  axis  of 
the  disk  is  given  in  argument  AR,  and  the  distance  z  from  the  plane  of  the  disk  is 
given  in  argument  AZ.  The  potential  of  the  disk  is  stored  in  function  FP. 

SUBROUTINE  CDSKO  (AA,  AR,  AZ,  FO) 

**#*******************************************************, ********************************** 
FORTRAN  SUBROUTINE  FOR  SOLID  ANGLE  OF  CIRCULAR  DISK 

**#*********4^**********#***#****#**#*#*¥*****4^********************************************* 


The  radius  a  of  the  disk  is  given  in  argument  AA,  the  distance  r  from  the  axis  of 
the  disk  is  given  in  argument  AR,  and  the  distance  z  from  the  plane  of  the  disk  is 
given  in  argument  AZ.  The  solid  angle  of  the  disk  is  stored  in  function  FO. 


RECTANGULAR  PLATE 


Analysis 


Let  2a  be  the  length  of  a  plate  and  let  26  be  the  breadth  of  the  plate.  Let  x,  y,  z 
be  Cartesian  coordinates  with  x  in  the  direction  of  the  length  of  the  plate,  with  y  in 
the  direction  of  the  breadth  of  the  plate,  and  with  z  in  the  direction  perpendicular 
to  the  plate.  Let  a  line  be  constructed  through  the  field  point  and  perpendicular  to 
the  plane  of  the  plate.  The  perpendicular  line  intersects  the  plane  of  the  plate  at  a 
point  with  coordinates  x,  y  with  respect  to  the  center  of  the  plate.  Let  u,  v  be  the 
coordinates  of  a  point  on  the  plate  with  respect  to  the  point  of  intersection.  Then  the 
potential  at  the  field  point  is  given  by  the  equation 


z) 


dudv 

Vu*  +V!  +  2* 


(54) 


Initial  integration  leads  to  an  inverse  hyperbolic  function,  then  final  integration  is 
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completed  with  an  integration  by  parts.  Introduction  of  limits  of  integration  leads  to 
the  equation 


ip  =  (a  -  x)sinh  =====  +  (a  -  x)sinh-1  -  ^—+  ^ 

v  (a  -  x)2  +  z2 

+  (a  +  x)sinh'*— =  ===  =  +  (a  +  x)sinh~‘ 


V(a  +  x)2  +  z2 

+  (6  -  y)sinh~‘  r7}a  XJ=  -  =  +  (b  -  yjsinh'1 


V(b  -  y)2  +  z2 


+  (b  +  y)sinh  1  -7^^  XJ-=  —  +  (6  +  y)sinh  *- 


V(a-x)2  +  z2 
( b  +  y) 

V(a  +  x)2  +  z2 
(a  +  x) 

^(6  -  y)2  +  z2 
(a  +  x) 


z  tan'1 


V(b  +  y)2  +  z2 
(a  -  x)(6  -  y) 


z  tan' 


z\/(a  -  x)2  +  (6  -  y)2  +  z2 
(a  +  x)(b  -  y) 


z  tan' 


z  tan' 


V(b  +  y)2  +  z2 

(a  -  x)(b  +  y)  _ 

zV(a  -  x)2  +  (b  +  y)2  +  z2 
(a  +  x)(b  +  y) 


, -  — - - —  -  <£■  ion  .  —  ■  -  ■  — ^ -  — 

zV (a  +  x)2  +  (b  -  y)2  +  z2  z\f ’ a  +  x)2  +  (b  +  y)2  +  z 


(55) 


Partial  differentiation  and  cancellation  of  terms  gives  the  components  of  the  gradient 
of  the  potential. 

Differentiation  with  respect  to  x  leads  to  the  equation 


- 4-sinh'1  ■  £.+?± 


dp  .  , 

—  =  sinh  - - _ - - 

dx  V(a  -  xTTP 


sinh 


(b  -  y) 


-  sinh 


(5  +  y) 


V(a  +  x)2  +  z2  V(o  +  x)2  +  z2 

differentiation  with  respect  to  y  leads  to  the  equation 


dip 

dy 


(a  -  x) 


(a  +  x) 


sinh 


V(b  -  y)2  +  z2 
.  (q  -  x) 


-  sinh  1 


V(b-y)2  +  z2 
(a  +  x) 


V(b  +  y)2  +  z2  V  (b  +  y)2  +  z2 

and  differentiation  with  respect  to  z  leads  to  the  equation 


(56) 


(57) 


dip  ,  (a-x)(b-y) 

—  =  tan  — r  — ■  -  -  -  — 1  ■  = 

dz  zV(q  -  x)2  +  (b  -  y)2  +  z2 


+  tan  1 


(a  -  x)(b  +  y) 


+  tan  1 


(a  +  x)(b  -  y) 


+  tan 


zV(a  -  x)2  +  (b  +  yj2  +  z2 
(a  +  x)(b  +  y) 


zV(a.  +  x)2  +  (b  -  y)2  +  z2  zv'(a  +  x)2  +  (b  +  y)2  +  z! 


(58) 


10 


Programming 

SUBROUTINE  RPLTP  (AA.  AB,  AX,  AY,  AZ,  FP) 

*#* ****************** ************ *********************** *******  *»«*«•****•**»«***** ********** 
FORTRAN  SUBROUTINE  FOR  POTENTIAL  OF  RECTANGULAR  PLATE 
********************************************************************************************* 

The  half-length  a  of  the  plate  is  given  in  argument  mA,  the  half- breadth  6  of  the 
plate  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  the  arguments  AX,  AY,  AZ.  The  potential  of  the  plate  is  stored  in  function  11 

SUBROUTINE  RPLTO  (AA,  AB.  AX,  AY.  AZ.  FO) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  SOLID  ANGLE  OF  RECTANGULAR  PLATE 
********************************************************************************************* 

The  half-length  a  of  the  plate  is  given  in  argument  AA,  the  half-breadth  6  of  the 
plate  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x.  y,  z  of  the  field  point 
are  given  in  the  arguments  AX,  AY,  AZ.  The  solid  angle  of  the  plate  is  stored  in  function  '  0. 


NONUNIFORM  PLATE 

Analysis 

The  potential  of  the  uniform  rectangular  plate  is  given  by  the  equation 


pb-v  _ dud 

J-b-  y  J -a-x  V  U®  +  V' 


+  Z 


and  the  derivatives  of  the  potential  are  given  by  the  equations 


dip 

p  +  b-y 

dx 

L-v 

p  +  b-v 

dy  J 

-b-v 

^  =  +  1 

p+b~v 

dz  J 

-b-v 

3 

,2)2 


.  dv 


3 


.  dv 


3 

,212 


(59) 

(60) 

(61) 

(62) 


If  the  surface  density  of  the  rectangular  plate  can  be  expressed  as  a  polynomial  in 
the  powers  of  *  +  u  and  y  +  v,  then  the  potential  of  the  nonuniform  plate  and  its 
derivatives  are  expressible  in  terms  of  members  of  a  family  of  integrals  of  which  the 
integral 


rw  r^^y^dudv 

-b -v  }U®  +  V*+Z*|* 


(63) 


is  the  m,  nth  member. 


1 1 


Integrations  of  lowest  degree  with  respect  to  u  are  given  by  the  equations 
du  u 


J 

J 

J 


ju2  +  t/2  +  Z2J^  (z*  +  V*)Vu*  +  V*  +  Z* 
(x  +  u)du 


XU 


f  U2  +  V*  +  zz]%  (2*  +  v2)Vul  +  vz  +  22  v'u2  +  V*  +  zl 

(x  +  u)zdu  xzu  2  x  +  u 


+  sinh  1 


j-u2  +  v2  +  Z2j^  (z2  +  H2)Vu2  +  tl2  +  Z2  V^U2  +  V*  +  2*  Vz'  +  U'1 


Vz 


(64) 

(65) 

(66) 


Integrations  of  higher  degree  are  given  by  the  recurrence  equation 
(x  +  u)mdu  (:r  +  u)m~l 


/ 


IL*  +  V*  +  Z 


2!2 


(m  -  2)Vtt2  +  v2  +  22 


(2m  -  3) 


(m  -  2) 


r  (x  +  ^ 

J  4-  ,, 


)m-ldu 


ju2  +  vz  +  z2| 2 
(r2  +  n2  +•  2 2) 


(m _ l)  („Z  J.  .,2  J.  ,2 

(m  -  2) 


f  +  u)m  *dli 

**  jli2  J-  V2  +  22!2 


Inasmuch  as  the  square  of  u  satisfies  the  identity 


v 2  +  22  -  i2  +  y2  +  22  -  2y(y  +  ti)  +  (y  *  v)2 


(67) 


(60) 


it  follows  that  the  recurrence  equation  replaces  the  integrals  of  lowest  degree  by  the 
products  of  power  polynomials  in  y  +  v  and  the  three  basic  functions 


u 


1 


(22  +  H2)\/u2  +  V‘  +  2‘  V  U‘  +■  V*  +  Z*  \  2*  ♦  tl* 


sinh 


(69) 


Thus  integration  with  respect  to  v  is  completed  with  the  aid  of  three  sets  of  integrals. 

The  integrations  of  lowest  degree  with  respect  to  v  for  the  first  set  are  given  by 
the  equations 


"J 

UI 


dv 


1 


(z2  +  n2)v/u2 

( y  +  v)dv 


tan' 


=  -  tan  '1- 


uv 


*  1  log 


v'u  +  v  z 


(22  +  t!2)  Vlt2  +  H2  +  Z2  2 
while  the  integrals  of  higher  degree  are  generated  by  the  recurrence  equation 
( y  +  v)ndv  f  (y  +  v)n~zdv 


(TO) 


(71) 


+  z‘ 


u  f  (y  +  V )ndv  -  U  f  +  V^r 

J  (z2  +  V2)V  U2  +  Vz  +  22  J  Vlt2  +  H 

f  (y  +  ti)nldn 

+  2yu  — - z. 

J  (22  +  ll2)v  U2  +  V 

-  (y2  +  z2)u 


+  Ze 


*J 


(y  +  v)n  2dv 


(z2  +  V2 


(72) 
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This  recurrence  is  cycled  in  the  direction  of  ascending  degree  if  the  arguments  satisfy 
the  stability  criterion 


(y2  +  z2)  S  3(y  +  t>)2 


(73) 


Otherwise  the  recurrence  is  cycled  in  the  direction  of  descending  degree  as  expressed 
by  the  equation 


u  1 

P  (y  +  v)ndv 

y2  +  z2  J 

1  \'u2  +  v2  +  z2 

2yu  | 

j"  (y  4  v)n*ldv 

y2  +  z2  J 

'  (z2  4  t;2)v/U2  4  V2  4  Z2 

f  (y  +  v)n*2dv 

y2  +  z2  J 

1  (Z2  4  V2)Vll2  4  V2  4  Z2 

(74) 


The  integrations  of  lowest  degree  with  respect  to  v  for  the  second  set  are  given  by 
the  equations 


J 

I 


dv 


'/u 


=  sinh  1  . — - — = 


+  v‘  +  z 
(y  +  v)dv 
v/u2 *■  v2  +  z2 


=  y  sinh"1- 


V  U2  +  Z2 


VU2  4  V2  4  Z2 


(75) 

(76) 


while  the  integrals  of  higher  degree  are  generated  by  the  recurrence  equation 


J 


{y  +  v)ndv  (y  +  v)n  1 


Vu2 


4  V  +  z 


2  n 

(2n-  1) 


v  u*  +  v  +  z‘ 


'dv 


+  z 


■  -  i)  r  ( y  +  v)n  'c 

n  J  Vtt2  +  V2  + 

(u2  +  y2  +  z2)  f  .<^3 

72  J  \  U2  V2  +  Z2 


(77) 


This  recurrence  is  cycled  in  the  direction  of  ascending  degree  if  the  arguments  satisfy 
the  stability  criterion 


(u2  +  y2  +  z2)  S  3(y  +  v)2 


(78) 


Otherwise  the  recurrence  is  cycled  in  the  direction  of  descending  degree  as  expressed 
by  the  equation 


J 


(y  +  v)ndv  (y  +  v)n*1 

Vu2  +  v2  +  z 2  (n  +  l)(u2  +  y2  +  z2) 


(2n 

(n  +  l)(u2 


(n  +  2) 
(n  +  l)(u2 


j*-  3)y _  r  (y+y^'dv 

2  4  y2  +  Z2)  J  V  U2  4  I/2  +  2 
!)  f  (y  +  v)n*2dv 

y2  +  z2)  J  \  u2  +  v2  +  z1 


(79) 
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The  integration  of  lowest  degree  with  respect  to  v  for  the  third  set  is  given  by  the 
equation 


f  sinh-1— 7=^ 

J  V? 


J  \fzz+vz  V 22  +  U2  J  Vuz+vz+zz  J  (z2+v2)Viz2+v2+z2 

while  the  integrals  of  higher  degree  are  generated  by  the  recurrence  equation 
f,  i  ,  u  (z2  +  v2)(  y  +  v)n~i  .  u 

J  <»  +  7P77* dv  ' - <^T) - °mh 


/  V? 


1  +  t/)"sinh  l- 


u  f  (y  +  v)ndv 

(n  +  1)  J  >/u2  +  v2  +  z2 

yu  r  (y  +  v )n~ldv 

(n.  +  1 )  J  vu!  +  v2  +  z2 


(n_ 


(n  +  1) 


( V 2  +  **) 


dv  (Bl) 


This  recurrence  is  cycled  in  the  direction  of  ascending  degree  if  the  arguments  satisfy 
the  stability  criterion 

(y2  +  z2)  §  3(y  +  v)2  (82) 

Otherwise  the  recurrence  is  cycled  in  the  direction  of  descending  degree  as  expressed 
by  the  equation 


(y  +  v)nsinh' 


_  .  _  (z2  +  v2)(y  +  v)n*‘  , 

~  o/v  .  . » <  a  2\  sinh 

^  (n  +  l)(y2  +  z2) 


_ yu _ 

(n  +  l)(y2  +  z2) 


(n  +  l)(y2  +  z2) 

2(n  +  2)y 
(n+  l)(y2  +  z2) 

(n  +  3) 

(n  +  l)(y2  +  z2) 


f  (y  +  v)n*ldv 

J  x/u^T V2  +  z2 

f  (y  +  v)n*zdv 

J  Vu2  +v2  +  z2 

j  (y  +  v)n+lsinh-1- 
J*  (y  +  u)n+2sinh_l- 


7=f=™j  ^  (83) 
'  z*  +-  v* 


The  cycling  of  each  recurrence  equation  in  the  direction  of  descending  degree  is 
started  at  the  64th  integral  in  order  to  achieve  full  accuracy  at  the  lowest  degree. 
Validity  of  all  formulae  for  indefinite  integrals  may  be  verified  directly  by  differentiation. 

Definite  integrals  are  evaluated  by  an  external  subroutine  which  refers  to  an  internal 
subroutine  for  the  evaluation  of  indefinite  integrals.  References  are  made  to  the 
internal  subroutine  with  arguments  set  equal  to  the  limits  of  integration  The  definite 
integrals  are  stored  in  a  matrix.  Accuracy  of  computation  can  be  verified  by  comparisons 


between  evaluations  by  subroutine  and  high-order  numerical  integrations,  if  the  value 
of  z  is  not  too  small. 

If  the  arguments  satisfy  the  inequality 

xz  +  yz  +  z2  ^  2 (a2  +  62)  (84) 


then  the  integrals  are  evaluated  with  16-point  Gaussian  integration.  The  m,  nth  member 
of  the  family  of  integrals, 


umvn  du  dv 
j(n  -  x)z  +  {v  -  y)z  + 


(85) 


is  derived  through  change  of  variable  from  x  +  u  to  u  and  from  y  +  v  to  v. 


Programming 

SUBROUTINE  RPLTM  (AA,  AB,  NA,  NB,  AX,  AY.  AZ,  FM) 

*********%*♦********************************************************************  +  ************ 
FORTRAN  SUBROUTINE  FOR  POWER  INTEGRATION  OVER  RECTANGULAR  PLATE 
********************************************************************************************** 

The  half-length  a  of  the  plate  is  given  in  argument  AA,  and  the  half-breadth  6  of 
the  plate  is  given  in  argument  AB.  The  number  of  powers  m  of  r+u  is  given  in  the 
argument  NA,  and  the  number  of  powers  n  of  y  +  v  is  given  in  the  argument  NB.  The 
Cartesian  coordinates  x,y,z  of  the  field  point  are  given  in  the  arguments  AX,  AY,  AZ. 
The  subroutine  constructs  three  arrays  of  integrals  with  respect  to  v,  then  synthesizes 
polynomials  to  complete  integration  with  respect  to  u.  The  definite  integrals  are  stored 
in  the  m  x  n  matrix  FM.  The  maximum  order  of  matrix  is  limited  to  32  '<  32. 


SPHERICAL  POLAR  COORDINATES 

Analysis 

Let  x.y.z  be  the  Cartesian  coordinates  of  a  point  in  space  and  let  i,j.  k  be  unit 
vectors  in  the  directions  of  increasing  r,  y,  z.  Let  r,  6,  if)  be  spherical  polar  coordinates 
and  let  ej,e2,e3  be  unit  vectors  in  the  directions  of  increasing  r,0,<f>  The  Cartesian 
coordinates  are  expressed  in  terms  of  the  spherical  polar  coordinates  by  the  equations 


x  r  sin  9  cos  0 

(86) 

y  -  r  sin  9  sin  <f> 

(87) 

z  -  r  cos  9 

(88) 

The  spherical  polar  coordinates  are  expressed  in  terms  of  the  Cartesian  coordinates 


15 


by  the  equations 


r  =  n/z2  +  yz  +  z2 

(89) 

.  Vxz  +  yz 

8  =  tan"1 - 

(90) 

z 

<p  =  tan"‘- 

(91) 

The  position  vector  r  is  given  by  the  equation 

r  =  r  sin  6  cos  <j>  i  +  r  sin  8  sin  0  j  +  r  cos  8  k  (92) 

and  the  differential  element  of  volume  is  given  by  the  equation 

|dri  =  r2  sin  8  dr  d.6  dip  (93) 

The  unit  vectors  e|,e2,e3  are  expressed  in  terms  of  the  unit  vectors  i.j.k  by  the 
equations 


e,  =  sin  8  cos  <p  i  t  sin  8  sin  <p  j  +  cos  8  k 

(9-4) 

e2  =  cos  8  cos  <p  i  *■  cos  8  sm  <p  j  -  sin  8  k 

(95) 

e3  =  -  sin  ip  i  +  cos  <P  ) 

(96) 

The  unit  vectors 

i,  j,  k  are  expressed  in  terms  of  the  unit  vectors  e, 

,  e2,e3  by  the 

equations 

i  =  sin  8  cos  0  e,  *■  cos  8  cos  <p  e2  sin  0  e3 

(97) 

j  =  sin  8  sin  0  ct  +  cos  8  sin  0  e2  +  cos  0  e3 

(98) 

k  =  cos  8  €  j  --  sin  8  ez 

(99) 

The  unit  vectors  are  right-handed  and  orthogonal.  Any  vector  or  tensor  invariant  can 
be  referred  interchangeably  to  either  set  of  unit  vectors.  The  components  of  the 
invariant  are  derived  from  the  scalar  products  of  the  invariant  with  the  unit  vectors 
Let  f  be  a  function  of  r,  8,  <p  as  expressed  by  the  equation 

ip  =  /(r,  8,  0)  (100) 

The  gradient  Vip  is  expressed  by  the  equation 


-  €i 


dip  e2  dip  e3  dip 
dr  r  90  r  sin  8  dip 


(101  ) 
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The  gradient  of  the  gradient  VVip  has  the  matrix 


VV0  = 


9zip 

1  dip\ 

1 

9f\ 

9r 

,r  90  ) 

9r\r 

sin  6 

d<t>) 

9 

( 1 

1 

dip 

i 

1 

9ip\ 

9r 

\  r  99/ 

r 

dr 

r2  99z 

99\r 

Jsin  0  d(f>  / 

1  dip\ 

9 

(- 

1  df\ 

1  dip  cos  9 

dip 

1  dzip 

9r\r 

sin  9  dip  ) 

99 

W2 

sin  0  dip  ) 

r  9r  r2sin  9  90 

r2sin20  902 

(102) 


The  trace  of  the  matrix  is  the  Laplacian, 

1 


V'V^  = 

Laplace’s  equation  is 


r2  dr 


(-5) 


a  j  .  av\ 

?  —  sin  5  — 

r2sm  0  90  \  90  J 


1  920 

r2sin20  9$ 2 


(103) 


VVf  =  0  (104) 

It  may  be  solved  by  the  method  of  separation  of  variables. 

Let  ip  be  represented  by  R®$  where  R  is  a  function  of  r  alone,  0  is  a  function  of  0 
alone,  and  4»  is  a  function  of  <f>  alone.  Substitution  in  Laplace’s  equation  shows  that 
the  factors  are  solutions  of  ordinary  differential  equations  which  are  linked  together 
through  arbitrary  constants  n,  m.  The  constants  n,  m  must  be  integers  in  order  that 
the  functions  shall  be  cyclical  with  respect  to  0,  <p. 

The  equation  for  R  is 


1  d  /  zdR\  n(n  +  1 ) 
ZR  dr\  dr) 


=  0 


The  function  R  is  any  linear  combination  of  the  functions 


1 

In+T 


The  first  derivatives  of  the  functions  are  given  by  the  equations 

d  (r")  =  nr"-‘ 


dr 


d  I  1  \  (n  + 


1) 


(105) 


(106) 


(107) 

(108) 


and  the  second  derivatives  of  the  functions  are  given  by  the  equations 

(r")  =  n(n-  1  )rn‘2 

arfc 

d2  /  1  (n+l)(n  +  2) 

dr2\rn+1  /  +  rnt3 


(109) 

(HO) 
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The  equation  for  8  is 


1  d  /  ■  n 

- sin  8  , 

sin  9  0  d8\  dB  / 


d@\ 


)  +  n(n  +  1 ) - — 


0 


sin20 

The  function  8  is  any  linear  combination  of  the  functions 

P™(cos  8)  Q?( cos  0) 


(HI) 


(112) 


where  P™(c os  0)  is  an  associated  Legendre  function  of  the  first  kind  and  Q™(cos  8)  is 
an  associated  Legendre  function  of  the  second  kind.  The  associated  functions  are 
defined  in  terms  of  the  regular  functions  by  the  equations 


Pn(co s  8)  -  sinm0 
Q™(cos  0)  -  sinm0 


An  analysis  of  Legendre  functions  is  given  in  Appendix  A. 

The  function  Pn(c os  0)  is  expressed  by  the  Rodrigues  formula 


dmPn(cos  0) 

(113) 

d(cos  0)m 

dmQn( cos  0) 

(114) 

d( cos  0)m 

Pn( cos  0) 


(  IT  dnsm2"0 
2nn!  d(cos0)n 


(115) 


The  function  Pn( cos  0)  is  a  power  polynomial  of  the  nth  degree  in  cos  0  The  function 
Qn( cos  0)  is  given  by  the  equation 


Qn( cos  0) 


|Pn(cOS0)log^~~-~)  Wn  ,  (COS  0) 


(116) 


where  Wn_,(cos  0)  is  a  polynomial  of  degree  n  1  in  cos  0.  The  functions  of  the  first 
kind  are  finite  everywhere  whereas  the  functions  of  the  second  kind  have  logarithmic 
singularities  at  the  poles.  When  sin  0  approaches  zero  the  associated  functions  of  the 
first  kind  approach  zero  like  sinm0  whereas  the  associated  functions  of  the  second 
kind  approach  infinity  like  cscm0.  When  sin  0  is  not  zero  both  functions  with  increasing 
n  approach  asymptotic  values  where  the  ratio  between  the  amplitude  of  the  second 
kind  and  the  amplitude  of  the  first  kind  is  |tt. 

The  functions  of  lowest  order  are  relatively  simple  and  the  functions  of  progressively 
higher  order  may  be  generated  from  them  with  the  aid  of  three- term  recurrence 
equations.  When  an  error  has  been  introduced  into  the  recurrence  at  the  Arth  cycle, 
it  may  be  represented  by  a  linear  combination  of  PJJ*  and  <?”  such  that  the  error  is 
equal  to  zero  for  the  ( k  l)th  cycle  but  is  equal  to  e  for  the  fcth  cycle.  Application 
of  the  recurrence  equations  to  the  linear  combination  of  PJJ*  and  QJJ*  changes 
progressively  the  order  of  each  term  If  the  recurrence  is  used  to  generate  P™,  then 
the  recurrence  must  be  cycled  in  that  direction  in  which  the  ratio  Q™/P™  decreases, 
or  if  the  recurrence  is  used  to  generate  <?£*,  then  the  recurrence  must  be  cycled  in 
that  direction  in  which  the  ratio  P™-  Q™  decreases.  Otherwise  the  relative  rounding 
error  will  not  remain  bounded. 


A 
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The  functions  of  lowest  order  are  given  by  the  equations 


,  /  1  •  cos  9  \ 

p.-i  «.  =  r.«.  ■  Si.g( O'?) 

=  cos  8  c,  =  P|Q0  -  1  (118) 

The  functions  for  m  =  0  then  satisfy  the  recurrence  equations 

nPn-i  -  (2 n  +  l)cos  0  Pn  +  (n  +  l)Pnn  ~  0  (119) 

n.Qn-l  -  (2n  +  l)cos  8  Qn  *-  (n  ?  1  )QnM  0  (120) 

Except  for  the  logarithmic  singularity  in  Qn  the  regular  functions  do  not  differ  greatly 
in  magnitude  and  both  recurrences  may  be  cycled  in  the  direction  of  increasing  n. 
The  associated  functions  satisfy  the  recurrence  equations 

(nf  m)P”,  -  (2n  +  t)cos0P™  -  (n  m-DP™,  0  (121) 

(n  +■  m)$™,  -  (2n  +  l)cos  8  Q™  f  (n  -  m  ~  1  )Q™, ,  0  (122) 

The  recurrence  becomes  more  sensitive  to  direction  with  increase  in  m  and  must  be 
cycled  toward  increasing  n  for  P™  but  toward  decreasing  n  for  Q™  In  the  ease  of  P ™ 
the  recurrence  is  started  with  the  equations 

P£.,(cos0)  =  O  (123) 

P£(cos  8)  =  sin^  (121) 

c  m: 


and  then  the  functions  are  generated  at  constant  m  for  progressively  increasing  n 
If  sin  9  is  very  small,  and  m  is  very  large  the  value  of  P™  can  be  below  the  index 
range  of  the  computer.  Such  a  situation  is  avoided  when  the  recurrence  is  applied  to 
the  derivatives  of  Pn  prior  to  multiplication  by  the  powers  of  sin  0. 

The  functions  with  common  n  satisfy  the  recurrence  equations 


(n  +  m)(n-m+  l)P?~l  -  2m  P?  -P™'1  0  (1251 

(n  +  m)(n-m+l)Q?-l-2m-~Q?  rQ?'1  0  (1261 

sin  8 


which  must  be  cycled  toward  decreasing  m  for  PJJ*  but  toward  increasing  m  for  Q”' 
In  the  case  of  <?£*  initial  values  are  derived  from  Q0.Qt,  through  the  use  of  the 
recurrence  equations 

rxQn  -i  _  (2n  +  1  )cos  fi  )  (n  ‘  1  )Qn. ,  0  (127) 

Q‘.l  +  (2n+l)sin0Qn-fl‘M  0  (1211) 

and  then  the  functions  are  generated  at  constant  n  for  progressively  increasing  m 
The  associated  function  of  the  first  kind  is  defined  by  the  equation 


(  1 29 ) 
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r 


The  first  derivative  of  the  associated  function  is  given  by  the  equation 

dP™  m  ,  dmPn 

—  =  m  sinm~‘0  cos  9  — - -r— 

dd  d( cos  0)m 


d(cos0)m+I 

and  the  second  derivative  of  the  associated  function  is  given  by  the  equation 

d2Pm  dmP 

—2  =  m(m  -  l)sin"-*fl  V 

d02  d(cos0)m 


(130) 


|m  -  n(n  +  l)jsinm0 


dmP„ 


d(co s  0)" 


sinm0  cos  9 


d( cos  0)T 


(131) 


Other  derivatives  which  are  important  for  the  computation  of  space  invariants  are 
given  by  the  equation 


de 


J  pm  \ 

(  — 2-  =  m(m  -  l)sinm  29  cos  9 
\sin  9  } 


dmP„ 


d(co s  9)v 


m  sinm0 


'P 

X  TV 


d(cos0)m+1 


(132) 


and  by  the  equation 


m2P™  cosB  dP ™  (  .  .  m  ,  dmP„ 

— -sr  +  — r  ™  m  -  l)sinm-20 

sm20  sin  9  d9 


d(cos  0)T 


-  m  sinm0 


dmP„ 


-  sinm0  cos  9 


d( cos  0)m 
dm*'P 


(133) 


d(cos  e)mM 

Wherever  negative  powers  of  sin  0  appear  in  the  formulae  they  are  multiplied  by  zero. 
The  equation  for  0  is 

1  d2<f> 

*d*2+m  ° 

The  function  0  is  any  linear  combination  of  the  functions 

cos  m0  sin  m<p 

The  first  derivatives  of  the  functions  are  given  by  the  equations 

d  cos  7710 
"  d0 

d  sin  m0 
d0 


=  -  m  sin  7Ti0 
-  +  m  COS  7710 


(13-1) 

(135) 

(136) 

(137) 
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and  the  second  derivatives  of  the  functions  are  given  by  the  equations 

cf2c os  7710  . 

- -z — ■  =  -  m  cos  mg> 

d02  Y 


d2sin  m<p 
ct02 

Let  Tp  be  expressed  by  the  equation 


V  . 


=  -  77l2sin  7710 


P™(co s  8)  eim * 


(138) 

(139) 


(140) 


All  77i  derivatives  of  Pn(cos  6)  are  computed  and  are  stored  in  an  array.  Initial 
summations  with  respect  to  n  establish  five  arrays  with  the  elements 


cm 

rn+l 

dmPn(<zos  8) 
cf(cos  0)m 

(141) 

nc™ 

rn+l 

dmPn( cos  0) 
d(cos  0)m 

(142) 

7i2cr 

dmPn{ cos  0) 

(143) 

rn  +  l 

d( cos  0)m 

Cm 

K'n 

rn+‘ 

dm+lPn(cos  0) 
d(cos0)m+1 

(144) 

< 

rn+l 

dm*'Pn{ cos  0) 
d(cos  0)m+1 

(145) 

Final  summations  with  respect  to  m  consist  of  complex  polynomial  evaluations  in 
powers  of  the  argument  sin  0  e x*. 


Programming 


SUBROUTINE  SPHPDV  (WO,  AR,  AQ,  AF,  NC,  SC,  CC,  RP,  PF.  DF,  DO) 

***************************************************************************************** 
FORTRAN  SUBROUTINE  FOR  SPHERICAL  POTENTIAL  AND  ANGULAR  DERIVATIVES 


********* 


The  mode  of  operation  is  given  by  WO.  The  radius  r  is  given  in  the  argument  AR,  the 
polar  angle  8  is  given  in  the  argument  AO,  and  the  azimuth  angle  0  is  given  in  the 
argument  AF.  The  order  of  the  matrices  of  coefficients  is  given  in  argument  \C,  the 
matrix  of  coefficients  for  sin  77i0  is  given  in  the  array  SC,  and  the  matrix  of  coefficients 
for  cos  7n0  is  given  in  the  array  CC.  The  matrix  of  Legendre  functions  is  stored  in  the 
array  pP.  In  each  matrix  the  rows  are  numbered  in  the  direction  of  increasing  n  and 
the  columns  are  numbered  in  the  direction  of  increasing  m  The  upper  right-hand 
half  of  each  matrix  is  padded  out  with  zeros.  The  potential  p  is  stored  in  the  function 
PF  if  the  mode  of  operation  WO  is  0.  The  potential  and  the  first  derivatives  dip  90, 
dip/dip  are  stored  in  the  function  DF  and  in  the  2-array  DF  when  the  mode  of  operation 
MO  is  1.  The  potential,  the  first  derivatives,  and  the  second  derivatives  d2<p  902,  d2yj  9090. 
ay/902  are  stored  in  the  function  PF,  in  the  2-array  DF,  and  in  the  3  array  DD  when 
the  mode  of  operation  MO  is  2. 
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SUBROUTINE  SPHPGD  (MO,  AR,  AO,  AF,  NC,  SC,  CC,  RP,  PF,  GF,  DF) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  SPHERICAL  POTENTIAL  AND  SPACE  GRADIENTS 
******  ********************************************************************** ******* ********** 

The  mode  of  operation  is  given  by  MO.  The  radius  r  is  given  in  the  argument  AR,  the 
polar  angle  0  is  given  in  the  argument  AQ,  and  the  azimuth  angle  <p  is  given  in  the 
argument  AF.  The  order  of  the  matrices  of  coefficients  is  given  in  the  argument  NC, 
the  matrix  of  coefficients  for  sin  m0  is  given  in  the  array  SC,  and  the  matrix  of 
coefficients  for  cos  mtp  is  given  in  the  array  CC.  The  matrix  of  Legendre  functions  is 
stored  in  the  array  RP.  In  each  matrix  the  rows  are  numbered  in  the  direction  of 
increasing  n  and  the  columns  are  numbered  in  the  direction  of  increasing  m.  The 
upper  right-hand  half  of  each  matrix  is  padded  out  with  zeros.  The  potential  <p  is 
stored  in  the  function  PF  if  the  mode  of  operation  MO  is  0,  The  potential  and  the 
gradient  V<p  are  stored  in  the  function  PF,  and  in  the  3-array  GF  when  the  mode  of 
operation  MO  is  1.  The  potential,  the  gradient,  and  the  gradient  of  the  gradient  VV^ 
are  stored  in  the  function  PF,  in  the  3-array  GF,  and  in  the  9-array  DF  when  the 
mode  of  operation  MO  is  2. 


ERROR  FUNCTION 


Analysis 

The  error  function  erf  z  is  defined  by  the  equation 

2  r  z 

erf  2  -  — e  u  du 
Vrr  Jo 

The  Dawson  integral  H(z)  is  defined  by  the  equation 

H(z) 


f 

Jo 


eu  du 


and  is  expressed  in  terms  of  the  error  function  by  the  equation 

Vrr 

H(z)  -  -  I--—  erf(+xz) 

The  conventional  Fresnel  integrals  C(v)  and  S(v)  are  defined  by  the  equation 


C(  v)  +  iS(v) 


r.* 

Jo 


du 


and  are  expressed  in  terms  of  the  error  function  by  the  equation 


C(v)  +  iS{v)  =  Lti  erf 


7TXlj 


(146) 


(147) 


(148) 


(149) 


(150) 


Expansion  of  the  exponential  function  in  series  and  term  by  term  integration  leads 
to  the  equation 


erf  z  =  —  r-  Y, 


vtt  rn-o  (2m  +  1  )m! 


(151) 
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which  expresses  the  error  function  as  an  absolutely  convergent  ascending  series.  The 
complex  Fresnel  integral  is  defined  by  the  equation 


E(z) 


(152) 


where  the  path  of  integration  lies  within  that  part  of  the  complex  plane  from  which 
the  positive  real  axis  is  excluded.  The  phase  of  z  is  limited  to  the  range  0  to  2rr,  and 
the  phase  of  z x/2  is  half  the  phase  of  z  There  are  convergent  series,  rational 
approximations,  and  asymptotic  series  for  the  complex  Fresnel  integral.  The  convergent 
series  is  given  by  the  equation 


^  m 

(2m  +  1  )m! 


(153) 


The  substitution  z  -»  -z2  converts  the  series  for  Fresnel  integral  into  the  series  for 
error  function  as  expressed  by  the  equation 

erf  z  =  1  -  iv'2£’(-z2)  (154) 

If  the  argument  x  +  iy  satisfies  the  inequality 

x2  +  y2  S  1  (155) 


or  both  of  the  inequalities 

1  <  x2  +  y2  <  38  x2  -  y2  +  0.256  x2y2  SO  (1 56) 

then  the  error  function  is  computed  with  the  ascending  series. 

If  the  argument  x  +  iy  satisfies  both  of  the  inequalities 

1  <  x2  +  y2  <  38  x2  ~  y2  *  0.256  x2y2  >0  (157) 


then  the  error  function  is  computed  with  the  rational  approximation  of  the  Fresnel 
integral.  The  error  function  is  expressed  by  the  equation 

-t2  18 

erf  z  =  1  -  2C .  ■  E  (158) 


f  <5t 


where  the  phase  of  z  is  limited  to  the  range  -|n  to  +|n,  and  the  positions  and  the 
residues  et  are  for  the  approximation  of  the  Fresnel  integral  by  sets  of  poles 
If  the  argument  x  +  iy  satisfies  the  inequality 

x2  +  y2  38  (159) 


then  the  error  function  is  computed  from  the  asymptotic  series.  Repeated  integration 
by  parts  leads  to  the  equation 


erf  z 


e  *2  "  (  l)m(2m)' 

^  p2m  i  Zm  *  1 

\  rr  m=o  ^  m  z 


(160) 


for  which  the  phase  of  z  is  limited  to  the  range  -  to  +^tt,  and  ,V  5  38. 

Extension  of  the  range  of  phase  beyond  these  limits  is  accomplished  with  the  aid 
of  the  equation 

erf(  z)  =  erf(z)  (161) 


\ 
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which  expresses  the  symmetry  of  the  error  function. 

Programming 

SUBROUTINE  CERF  (MO,  AZ.  EF) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  COMPLEX  ERROR  FUNCTION 
********************************************************************************************* 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  argument  z 
are  given  in  array  AZ.  The  complex  error  function  is  computed  by  series  expansions 
and  rational  approximations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function 
erf  z  are  stored  in  array  EF.  If  MO  =  1.  and  the  phase  of  z  is  in  the  range  to  +|tt, 
the  real  and  imaginary  parts  of  the  function  1  -  erf  z  are  stored  in  array  EF. 


COMPLEX  GAMMA  FUNCTION 


Analysis 

The  gamma  function  F(z)  is  defined  by  the  equation 


* 


71  J 


(162) 


where  y  is  Euler’s  constant.  The  gamma  function  has  poles  at  the  negative  integers 
such  that  the  residue  of  the  nth  pole  is  (-1  )n/n!. 

For  a  small  argument  the  reciprocal  of  the  gamma  function  is  given  by  the  Bourguet 
convergent  series  and  for  a  large  argument  the  logarithm  of  the  gamma  function  is 
given  by  the  Stirling  asymptotic  series.  Intermediate  regions  can  be  spanned  by 
recurrence  relations.  A  rational  approximation  is  not  necessary. 

The  gamma  function  of  an  argument  with  a  negative  real  part  is  expressed  in  terms 
of  the  gamma  function  of  an  argument  with  a  positive  real  part  by  the  reciprocal 
equation 


r(z)r(i  -  z)  = 


7T 

sin  ttz 


(163) 


It  is  necessary  to  evaluate  series  expansions  only  for  arguments  with  positive  real 
parts. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yzSl  (164) 


then  the  gamma  function  is  derived  from  an  ascending  power  series.  The  reciprocal 
of  the  gamma  function  is  given  by  the  equation 


1 

r(i  +  z) 


00 

=  E  cmzm 


(165) 


for  which  the  coefficients  cm  are  derived  in  Appendix  B. 
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If  the  argument  x  +  iy  satisfies  the  inequality 

i2 +  yz g 32  (166) 

then  the  gamma  function  is  computed  from  a  descending  power  series.  From  the 
equations  on  page  252  of  reference  1,  the  logarithm  of  the  gamma  function  is  given 
by  the  equation 

log  r(2)  M*  -  i)  log  *  -  *  +  I  log  (2n)  4  (,6V) 

for  which  the  Bernoulli  numbers  BZm  are  derived  in  Appendix  B.  Summation  of  the 
series  is  continued  until  there  is  no  change  in  sum  or  until  m  ~  18. 

If  the  argument  x  +  iy  satisfies  the  inequality 

1  <  x*  4  yz  <  32  ( 1 68) 

then  the  gamma  function  is  computed  with  the  aid  of  the  difference  equation 

r(l4z)  =  zr(z)  (169) 

If  n  is  the  integer  which  is  nearest  in  value  to  x  and  if  n  satisfies  the  inequality 

)z  -  n!a  g  1  (170) 


then  the  gamma  function  is  given  by  the  equation 

F(z)  =  (z  -  !)•  "(z  -  n  4  l)r(z  -  n  4  1) 


(171) 


for  which  F(z-n4  1)  is  evaluated  from  the  convergent  series.  If  n  is  the  smallest 
integer  which  satisfies  the  inequality 


'  z  ^  n'z  S  32 

then  the  gamma  function  is  given  by  the  equation 

r(z  4  n) 

(z  4  n  -  1) 


r <*>  -  ^ 


(172) 


(173) 


for  which  T(z  4  n)  is  evaluated  from  the  asymptotic  series. 

Programming 

SUBROUTINE  CGAWMA  (VO,  AZ,  :  G) 

**#***♦********+***************»***  +  **********  +  **'***♦  +  *»**•#■*****♦*#*****************»***»**** 
FORTRAN  SUBROUTINE  FOR  COMPLEX  GAMMA  FUNCTION 

*******************************  *<r****-**-******:**************************************t,***.****»* 


The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  complex 
argument  z  are  given  in  array  AZ.  The  complex  gamma  function  is  computed  by  series 
expansions  and  recurrence  relations.  If  MO  -  0,  the  real  and  imaginary  parts  of  the 
complex  function  F(z)  are  stored  in  array  CC.  If  VO  =  1,  the  real  and  imaginary  parts 
of  the  complex  function  logF(z)  are  stored  in  array  rG. 

SUBR0UT,Nr  OGAVVA  (VO,  A"  FG) 

*+***********♦*****************************************************************************♦* 

FORTRAN  SUBROUTINE  FOR  DOUBLE-  PRECISION  GAMMA  FUNCTION 
********************************************************************************************* 


The  mode  of  operation  is  given  in  VO.  The  real  and  imaginary  parts  of  the 
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double -precision  argument  *  are  given  in  array  AZ.  The  complex  gamma  function  is 
computed  by  series  expansions  and  recurrence  relations.  If  WO  =  0,  the  real  and 
imaginary  parts  of  the  double-precision  function  r(z)  are  stored  in  array  FG.  If  WO  =  1, 
the  real  and  imaginary  parts  of  the  double-precision  function  log  r(a)  are  stored  in 
array  FG. 


COMPLEX  DIGAMMA  FUNCTION 


Analysis 

The  digamma  function  +(2)  is  defined  in  terms  of  the  gamma  function  T(2)  and  the 
derivative  of  the  gamma  function  T\z)  by  the  equation 


+(z) 


r(z) 

r(z) 


(174) 


For  a  small  argument  the  reciprocal  of  the  gamma  function  is  given  by  the  Bourguet 
convergent  series  and  for  a  large  argument  the  logarithm  of  the  gamma  function  is 
given  by  the  Stirling  asymptotic  series.  The  derivative  of  the  gamma  function  is  derived 
by  differentiation  of  the  series 

The  digamma  function  of  an  argument  with  a  negative  real  part  is  expressed  in 
terms  of  the  digamma  function  of  an  argument  with  a  positive  real  part  by  the 
reciprocal  equation 

+(2)  =  *(1  -  z)  -  rrcot  rrz  (175) 


It  is  necessary  to  evaluate  series  expansions  only  for  arguments  with  positive  real 
parts. 

If  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  yz%  1  (176) 

then  the  digamma  function  is  given  by  the  equation 


'K  1  +  z)  =  -  T(  1  +2)  £  77lCm2m“I 


for  which  the  coefficients  cm  are  derived  in  Appendix  B. 
If  the  argument  2  +  iy  satisfies  the  inequality 

xz  +  yzi32 

then  the  digamma  function  is  given  by  the  equation 

1  *  B 


'l'(z)  =  log  2 


X '  ^2m 

t,  2mzzm 


(177) 


(178) 


(179) 


for  which  the  Bernoulli  numbers  BZm  are  derived  in  Appendix  B.  Summation  of  the 
series  is  continued  until  there  is  no  change  in  sum  or  until  m=  IB. 

If  the  argument  x  +  iy  satisfies  the  inequality 

1  <.  xz  +  yz  <  32  (180) 


then  the  digamma  function  is  computed  with  the  aid  of  a  difference  equation.  If  n  is 
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(161) 


the  integer  which  is  nearest  in  value  to  x  and  if  n  satisfies  the  inequality 

| z  -  n|2  §  1 

then  the  digamma  function  is  given  by  the  equation 

*(*)  =  -—  +  -  + - - - : r+*(«-  n  +  1)  (182) 

2-1  2— n+1 

for  which  +(2-71+  1)  is  evaluated  from  the  convergent  series.  If  n  is  the  smallest 
integer  which  satisfies  the  inequality 

(z  +  n|2  £  32  (183) 

then  the  digamma  function  is  given  by  the  equation 

♦(«)=- i  -  — - - - + 'J '(z+n)  (184) 

2  2  +  n  -  1 

for  which  +(2  +  n)  is  evaluated  from  the  asymptotic  series. 

Programming 

SUBROUTINE  CPSI  (MO,  AZ,  PS) 

***********  **************************************************** *********************** ******* 
FORTRAN  SUBROUTINE  FOR  COMPLEX  DIGAMMA  FUNCTION 
********************************************************************************************* 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
2  are  given  in  array  AZ.  The  complex  digamma  function  is  computed  by  series  expansions 
and  recurrence  relations.  Calls  are  made  when  necessary  to  Subroutine  CGAMMA.  If 
MO  =  0,  the  real  and  imaginary  parts  of  the  function  F'(2)  are  stored  in  array  PS.  If 
MO  =  1,  the  real  and  imaginary  parts  of  the  function  +(2)  are  stored  in  array  PS. 


FIRST  ORDINARY  BESSEL  FUNCTION 


Analysis 

The  ordinary  Bessel  function  of  the  first  kind  Jn(z)  is  given  by  the  absolutely 
convergent  series  in  the  equation 


J»{*)  -  Z 


ml  l ^\n+2m 


(-0m(b) 


(185) 


l!0  m!(n  +  m)! 

The  convergent  series  is  used  if  the  order  and  the  argument  satisfy  the  criterion 

M  i  j|2)2  (186) 


A  descending  recurrence  is  used  to  extend  the  range  of  orders  to  lower  orders. 
The  Bessel  function  is  given  by  the  equation 


Jn(z)  =  ~  je2 nWiKn(iz )  -  e^nTr</fn(-i2) 


(187) 


in  which  the  modified  Bessel  functions  of  the  second  kind  can  be  expressed  by  rational 
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approximations.  Thus  the  Bessel  function  is  given  by  the  equation 

Jn(z)  =  — - — i  jl  +  E  — t  i|cos(z  -  |rvn  -  Jw)  +  isin(z  -  |nn  -  ~7t)| 

(2nz)*  '  »-»-«-«*)(  ) 

+  — — — i  j  1  +  E  - —  ]  )cos(z  -  |n7T  -  ^7t)  -  i  sin(z  -  knn  -  |tt)|  (180) 

(2nz)2  (  *=>  +  12  -  <5*  )  (  ) 


where  the  phase  of  z  is  limited  to  the  range  -|jt  to  +|tt,  and  the  positions  6k  and  the 
residues  et  are  for  the  approximation  of  the  modified  Bessel  function  by  sets  of  poles. 
The  rational  approximation  is  used  if  the  argument  x  +  iy  satisfies  the  inequalities 
i 

(x2  +  y2)*  §  17.5  -  |y|  +  0.096  x2  >  0  (189) 


The  rational  approximation  is  available  for  orders  0  and  1. 
The  Bessel  function  is  given  by  the  equation 


•Ai(z)  =  (~)  j^n(z)cos(z  -  |nrr  -  |n)  -  Qn(z)sin(z  -  ^nn  ~  |jr)| 


(190) 


where  Pn(z)  is  the  sum  of  the  even -ordered  terms  and  Qn(z )  is  the  sum  of  the 
odd-ordered  terms  in  the  asymptotic  series  in  the  equation 


P  (2)  ►  ,«  (2)  =  V  1^1:  (^  (^(m 


I)2) 


(191) 


with  N  §  36.  The  asymptotic  series  is  used  if  the  order  and  the  argument  satisfy  the 
criterion 


|z|  17.5  +  |n2  (192) 

An  ascending  recurrence  is  used  to  extend  the  range  of  order  to  larger  orders  if  the 
order  satisfies  the  criteria 


|z|  i  17.5  (193) 

and 

v'2(]z!  -  17.5)  §  |nl  §  ||z|  -  ||,7m  z!  +  |!||z|  -  [ 9m  z||  (194) 

Otherwise  the  ascending  series  is  used  with  a  descending  recurrence. 

Use  of  the  series  expansions  is  limited  to  positive  orders  and  to  arguments  with 
phase  in  the  range  -|tt  to  +|tt.  Extension  of  ranges  of  order  and  phase  is  accomplished 
with  the  aid  of  the  equations 

U.n(z)  =  (-  l)Vn(z)  =  Jn(- z)  (195) 

which  express  the  symmetry  of  the  Bessel  function. 

Programming 

SUBROUTINE  BSSLJ  (AZ.  IN,  FJ) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ORDINARY  BESSEL  FUNCTION  OF  INTEGRAL  ORDER 
********************************************************************  *****************  ******** 


The  real  and  imaginary  parts  of  the  argument  z  are  given  in  array  AZ,  and  the 


integer  order  n  is  given  in  IN.  The  ordinary  Bessel  function  of  the  first  kind  is  computed 
by  series  expansions  and  recurrence  relations.  The  real  and  imaginary  parts  of  the 
function  Jn(z)  are  stored  in  array  FJ. 


SECOND  ORDINARY  BESSEL  FUNCTION 


Analysis 

The  ordinary  Bessel  function  of  the  second  kind  Kn(z)  is  given  by  the  absolutely 
convergent  series  in  the  equation 

r„(*)  =  -  -1  V  (n~m,~1)!  (h)-nt2m 


71  TTX^O 


l  E  U 

71  m-0  L 


>nl  1 


y  +  log(|z)  -  -  E  -k  2  E  i 

m-0  L  ^  k=l  *  c  K 


(-i)">(§z)"*2 


(§z)n*Zm 

m!(n  -*  m)1 

The  convergent  series  is  used  if  the  argument  x  +  iy  satisfies  the  inequality 

l 

or  both  of  the  inequalities 

x*  +  y*<  289  -  \y\  +  0.096  xz  S  0 


(196) 


(197) 


(198) 


The  evaluation  of  the  convergent  series  for  Fn(z)  is  continued  to  convergence  of  the 
associated  series  for  Jn(z). 

The  Bessel  functions  are  given  by  the  equations 


Jn(*)  = 

i 

+  — 

77 

pn'Vn(iz)-  e  ^nniKn(-  iz)| 

(199) 

>«(*)  = 

1 

7T 

je^'V^iz)  +  e‘ z"’,',A'n(-tz)| 

(200) 

in  which  the  modified  Bessel  functions  of  the  second  kind  can  be  expressed  by  rational 
approximations.  Thus  the  Bessel  functions  are  given  by  the  equations 


1  (  14 

-  + - 1  )  1  +  E 

(2ttz)=  ' 

it 


iz  ~  <5* 


(cos(z  -  |n7T  -  \n)  i  sin(z  -  |nn  ‘n)| 


(201) 


(2rrz) 


(2tt  z) 


- — r  j  1  +  E  — r~*  — —  !  I  cos(z  -  inn  -  in)  -  i  sin(z  -  inn  -  in 
T~)i  (  +tz-6t)( 

SI 

1 1  cos (z  -  \nn  -  -7t)  -  i  sin(z  -  \m r  -  jTt)|  (202) 


(2tts) 

>'n(z)  = - ~  i  \  1  +  E  |  lcos(z  \nn  -  Jtt)  +  i  sin(z  --  \m r  - 

\  k=l  lZ  '  *k  )  (  ) 

14 

r  ^  E  -- 


iz  -  6t 


where  the  phase  of  z  is  limited  to  the  range  |tt  to  and  the  positions  6k  and  the 
residues  e*  are  for  the  approximation  of  the  modified  Bessel  function  by  sets  of  poles 
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The  rational  approximation  is  used  if  the  argument  x  +  iy  satisfies  the  inequalities 
xz  +  yz  <  289  -  |y|  +  0.096  xz  >  0  (203) 

The  rational  approximation  is  available  for  orders  0  and  1. 

The  Bessel  function  is  given  by  the  equation 
i 

Yn(z)  =  jpn(z)sin(z  -  -  |tt)  +  Q„(z)cos(z  -  \ rnr  -  Jn)|  (204) 

where  Pn(z)  is  the  sum  of  the  even-ordered  terms  and  Qn(z)  is  the  sum  of  the 
odd-ordered  terms  in  the  asymptotic  series  in  the  equation 


Pjz)  +  iQn(z)=  £ 

m=0 


ln*-(J)*}-|n»-(m-i)*| 


(205) 


m!(-2i2)m 

with  N  §  36.  The  asymptotic  series  is  used  if  the  argument  x  *  iy  satisfies  the  criterion 

xz  +  yz  Z  289  (206) 


The  rational  approximation  and  the  asymptotic  approximation  are  used  only  for 
arguments  with  phase  in  the  range  -|w  to  +^n.  Extension  of  phase  to  other  ranges 
is  accomplished  with  the  aid  of  the  equations 

Y„(z)  =  e'nvi  Yn(ze~ni)  +  2i  cos  nit  Jn(ze-ni)  (207) 


Yn(z)  =  etn,r‘  Yn(ze+ni)  -  2i  cos  nrr  J^ze**')  (208) 

where  the  first  equation  is  used  if  2  is  in  the  second  quadrant  and  the  second  equation 
is  used  if  2  is  in  the  third  quadrant. 

The  series  evaluations  are  used  only  for  orders  zero  and  one.  and  an  ascending 
recurrence  is  used  if  the  order  is  greater  than  one.  The  extension  of  order  to  negative 
orders  is  achieved  with  the  aid  of  the  equation 

Y_n(z)  -  (~l)nYn(z)  (209) 

which  expresses  the  symmetry  of  the  Bessel  function. 

Programming 

SUBROUTINE  BSSL.r  (AZ,  IN,  rY) 

********************************************************************************************* 
FORTRAN  SUBROUTINE  FOR  ORDINARY  BESSEL  FUNCTION  OF  INTEGRAL  ORDER 
********************************************************************************************* 

The  real  and  imaginary  parts  of  the  argument  z  are  given  in  array  AZ,  and  the 
integer  order  n  is  given  in  IN.  The  ordinary  Bessel  function  of  the  second  kind  is 
computed  by  series  expansions,  rational  approximations,  and  recurrence  relations. 
The  real  and  imaginary  parts  of  the  function  t'n(2)  are  stored  in  array  r  -. 


FIRST  MODIFIED  BESSEL  FUNCTION 


Analysis 

The  modified  Bessel  function  of  the  first  kind  In(z)  is  given  by  the  absolutely 
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convergent  series  in  the  equation 


In(z)  =  E 


”  (|z)"+2m 


(210) 


m=0  m!(n  +  m)! 

The  convergent  series  is  used  if  the  order  and  the  argument  satisfy  the  criterion 

|niSi|2l2  (211 

A  descending  recurrence  is  used  to  extend  the  range  of  order  to  lower  orders. 

The  Bessel  function  is  given  by  the  equations 


/,(*)  =  +  -  Wjze™)  -  e-^Kn(z) 


In(z) 


"  - ;  K 


(ze~ni)  -  enwiKn(z) 


(212) 

(213) 


where  the  first  equation  is  used  if  z  is  in  the  fourth  quadrant  and  the  second  equation 
is  used  if  z  is  in  the  first  quadrant.  The  modified  Bessel  functions  of  the  second  kind 
can  be  expressed  by  rational  approximations.  Thus  the  Bessel  function  is  given  by  the 
equation 


/„(*)  = 


z  -  6k 


-z±rvni±~n\  _ 

e  2  f 

+ - .  -•  <\ 


14 

V 


(27T2)2 


z  -  6t 


(214) 


where  the  ±  sign  is  the  same  as  the  sign  of  y.  and  the  positions  6k  and  the  residues 
ek  are  for  the  approximation  of  the  modified  Bessel  function  by  sets  of  poles.  The 
rational  approximation  is  used  if  the  argument  x  +  iy  satisfies  the  inequalities 


(xz  +  yz)z  S  17.5 


-  \x,  +  0.096v2  >  0 


The  rational  approximation  is  available  for  orders  0  and  1. 

The  Bessel  function  is  given  by  the  asymptotic  series  in  the  equation 

In{2)  a  _£l.  v‘  irt2  -(!)*!  •  in2  -(m  -  |)2| 

(27T2)2  m'° 


- -  y 

(2ttz)2  m=° 


m<(~2z)m 

(m-  |)2i 


m\{2z)T 


(215) 


(216) 


for  which  .V  S  36.  The  asymptotic  series  is  used  if  the  order  and  the  argument  satisfy 
the  criterion 


2^  17.5  +  in2  (217) 

An  ascending  recurrence  is  used  to  extend  the  range  of  order  to  larger  orders  if  the 
order  satisfies  the  criteria 


and 


iXe  z |  S  17.5 


(218) 


v  2(12]  17.5)  ^  |n!  ■-  ||z|  -  || 3te  zl  4  IXe  z\ 


(219) 
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V 


•v 


t  1 


Otherwise  the  ascending  series  is  used  with  a  descending  recurrence. 

Use  of  the  series  expansions  is  limited  to  positive  orders  and  to  arguments  with 
phase  in  the  range  -|tt  to  +|tt.  Extension  of  ranges  of  order  and  phase  is  accomplished 
with  the  aid  of  the  equations 

/_„(2)M-l)"/n(-2)  =  /n(z)  (220) 

which  express  the  symmetry  of  the  Bessel  function. 

Programming 

SUBROUTINE  BSSLI  (MO,  AZ.  IN,  FI) 

***************************************************** **************************************** 
FORTRAN  SUBROUTINE  FOR  MODIFIED  BESSEL  FUNCTION  OF  INTEGRAL  ORDER 
********************************************************************************************* 

The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  argument  2 
are  given  in  array  AZ,  and  the  integer  order  n  is  given  in  IN.  The  modified  Bessel 
function  of  the  first  kind  is  computed  by  series  expansions  and  recurrence  relations. 
If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  In{z)  are  stored  in  array  FI.  If 
MO  =  1,  and  the  phase  of  2  is  in  the  range  -  |rr  to  +^tt,  the  real  and  imaginary  parts 
of  the  function  e~*In(z)  are  stored  in  array  FI. 


SECOND  MODIFIED  BESSEL  FUNCTION 

Analysis 

The  modified  Bessel  function  of  the  second  kind  Kn(z)  is  given  by  the  absolutely 
convergent  series  in  the  equation 

n- 1 


*«(*)  =  3  2  ■  -1)-(7L,m  1)! 


-M)n  £ 


m! 

.  t  .  1  £  1  1  n"m  1 

J  +  !og(|z)  -  £  ~ 

C  j  fv  C  £ _  j 


(|2)n^m 

m!(n  +  m)! 


The  convergent  series  is  used  if  the  argument  x  +  iy  satisfies  the  inequality 

xz  +  y*  5  1 

or  both  of  the  inequalities 

xz  +  yz  <  289  +  x  +  0.096  yz  S  0 


(221) 


(222) 


(223) 


The  evaluation  of  the  convergent  series  for  Kn{z )  is  continued  to  convergence  of  the 
associated  series  for  In(z). 

The  Bessel  function  is  given  by  the  rational  approximation 


(224) 


for  which  the  positions  6t  and  the  residues  are  for  the  approximation  of  the  modified 
Bessel  function  by  sets  of  poles.  The  rational  approximation  is  used  if  the  argument 


32 


(225) 


x  +  iy  satisfies  the  inequalities 


xz  +  yz  <  289 


+  x  i-  0.096  yz  >  0 


The  rational  approximation  is  available  for  orders  0  and  1. 

The  Bessel  function  is  given  by  the  asymptotic  series  in  the  equation 


i 


n(  ’  \2z)  ho  m\(2zr 


I)21 


(226) 


for  which  N  ^  37.  The  asymptotic  series  is  used  when  the  argument  x  +  iy  satisfies  the 
inequality 


xi  +  yzi289  (227) 

The  series  evaluations  are  used  only  for  orders  zero  and  one,  and  an  ascending 
recurrence  is  used  if  the  order  is  greater  than  one.  The  extension  of  order  to  negative 
orders  is  achieved  with  the  aid  of  the  equation 

K.  n(z)  =  Kn(z)  (228) 

which  is  an  identity  for  all  orders  integer  or  complex. 


Programming 

SUBROUTINE  bsslk  (MO,  AZ,  IN,  FK) 

FORTRAN  SUBROUTINE  FOR  MODIFIED  BESSEL  FUNCTION  OF  INTEGRAL  ORDER 


The  mode  of  operation  is  given  in  MO.  The  real  and  imaginary  parts  of  the  argument 
z  are  given  in  array  AZ,  and  the  integer  order  n  is  given  in  IN.  The  modified  Bessel 
function  of  the  second  kind  is  computed  by  series  expansion,  rational  approximation, 
and  recurrence  relations.  If  MO  =  0,  the  real  and  imaginary  parts  of  the  function  A'n(z) 
are  stored  in  array  rK.  If  MO  1,  the  real  and  imaginary  parts  of  the  function  e*An(z) 
are  stored  in  array  FK. 


COMPLEX  BESSEL  FUNCTION 


Analysis 

Bessel  functions  of  complex  order  v  and  complex  argument  z  are  expressed  by 
absolutely  convergent  series.  The  Bessel  function  Jv(z )  is  given  by  the  equation 


.  ,  ,  _  y  (  ;l)m(^r2m 

v  ’  rto  m!r(u  +  m  f  1) 
and  the  Weber  function  t\,(z)  is  given  by  the  equation 

Jv{z)  cos  urt  -  J_v{z) 


= 


Sin  UTT 


(229) 


(230) 


If  v  is  a  negative  integer  n  then  all  terms  with  m<n  are  zero  in  the  series  for  ^.n(z) 
because  the  gamma  function  of  a  negative  integer  is  infinite.  Thus  the  Bessel  functions 


33 


satisfy  the  equation 


(231) 


J-niz)  =  (-l)n-Ai(z) 

and  the  Weber  function  is  given  by  the  equation 

Yn(z)  =  iim  Yv(z)  (232) 

l fix 

Tiie  Weber  function  Yn(z)  must  be  expressed  by  a  special  series. 

If  z  is  replaced  by  e±J,xz  in  the  convergent  series  then  Jv(z)  is  replaced  by  exv"xJ  v(z) . 
The  phase  of  z  is  limited  to  the  range  from  -it  to  +n  in  the  evaluation  of  J„(z),  but 
the  factors  et',wx  may  be  applied  to  Jt,( z)  in  order  to  extend  the  range  of  the  phase 
of  z  outside  the  range  from  -n  to  t7t.  If  v  is  real  the  absolute  magnitude  of  etvx,x  is 
unity,  but  if  v  has  any  imaginary  part,  then  the  magnitude  of  e*1'*'  may  be  small  or 
large  according  to  the  sign  of  the  imaginary  part  of  v. 

The  ratio  between  the  absolute  value  of  the  771th  term  and  the  absolute  value  of 
the  (m  -  l)th  term  is  given  by  the  expression 


m\v  +  m, 


(233) 


The  ratio  is  unity  wherever  the  terms  in  the  series  have  a  minimum  or  a  maximum. 

If  v  is  negative  and  real,  the  terms  of  the  convergent  series  may  increase,  decrease, 
increase,  and  decrease  with  increasing  order  r.i.  The  value  of  m  for  a  unit  ratio  between 
terms  is  estimated  by  the  equation 


(m  -  u  )  (234) 


and  by  the  equation 


rn 


2 


(m  ■  v  )  (235) 


The  requirement  that  m  be  real  and  positive  limits  the  number  of  minima  to  one  and 
the  number  of  maxima  to  two  unless  12  '•  \u:. 

If  v  is  positive  and  real,  the  terms  of  the  convergent  series  may  increase  before 
they  finally  decrease  with  increasing  order  m.  The  value  of  rn  for  a  unit  ratio  between 
terms  is  established  by  the  equation 

''.Pi  t  v  <p)*  -  2  2 

rn  - - -  (236) 


The  requirement  that  m  be  positive  limits  the  number  of  maxima  to  one  When  v  is 
complex  the  values  of  m  become  the  roots  of  a  quartic  equation. 

The  absolute  rounding  error  in  the  convergent  series  is  determined  by  the  relative 
rounding  error  in  the  largest  term  of  the  series.  In  order  to  control  the  relative 
rounding  error  in  the  series  itself  it  is  necessary  to  limit  the  evaluation  of  the  scries 
to  conditions  where  the  terms  diminish  after  the  first  term.  The  full  accuracy  of  the 
computer  is  guaranteed  only  if  the  order  and  the  argument  satisfy  the  criterion 

1  „  z 


A  practical  accuracy  of  computation  still  may  be  achieved  if  this  criterion  is  relaxed. 


especially  when  the  order  is  negative  or  the  argument  is  imaginary. 

The  Bessel  function  is  given  by  the  equation 
1 

Jv(z)  =  (~^\  |^i/(z)cos(z  -  |l/7!  -  in)  -  $„(z)sin(z  -  |t /n  -  ~n)|  (238) 


where  the  function  Pv( z)  is  the  sum  of  the  even-ordered  terms  and  the  function  <?„(z) 
is  the  sum  of  the  odd -ordered  terms  in  the  asymptotic  series 


Pv(z)  +  iQ„(z)  = 


V  n*  +  m  +  k) _ 

mir(i/-m  + |)(-2tz)m 


(239) 


The  ratio  of  gamma  functions  is  computed  from  the  product 


T(i/  +  m  + 
f(w  -  m  + 


i  \  m 

I-  =nwz-(k 


m 


(240; 


If  v  is  half  an  odd  integer  the  series  terminates  after  a  finite  number  of  terms. 

If  z  is  replaced  by  etZnlz  in  the  asymptotic  series,  then  the  terms  of  the  series  are 
reversed  in  sign,  whereas  the  actual  value  of  d„(z)  may  be  smaller  or  larger  according 
to  the  sign  of  the  imaginary  part  of  v.  This  failure  of  the  asymptotic  series  is  related 
to  the  Stokes  phenomenon.  Its  effect  is  diminished  if  use  of  the  asymptotic  series  is 
limited  to  arguments  with  phase  in  the  range  -g7T  to  +  |n,  whence  the  asymptotic 
series  is  corrected  by  the  factor  e11"1  when  the  actual  phase  is  outside  the  range 
from  -gtt  to  +  §7T. 

If  v  is  real  the  terms  of  the  asymptotic  series  may  increase,  decrease,  and  increase 
with  increasing  order  to.  The  ratio  between  the  absolute  value  of  the  mth  term  and 
the  absolute  value  of  the  (m  l)th  term  is  given  by  the  expression 


m’|2z! 


(241) 


The  ratio  is  unity  wherever  the  terms  in  the  series  have  a  minimum  or  a  maximum. 
When  v  is  .real,  the  value  of  m  for  a  unit  ratio  between  terms  is  estimated  by  the 
equation 

to  -  |  -  -  z;  ±  \  zi2  :z:  +  V®  (to  -  j  5?  v  )  (242) 

and  by  the  equation 

TO  \  =  +  ’z:  ±  \  i Z|®  *  z  +  |l/l*  (to  -  I  £  ,v')  (243) 


The  requirement  that  to  be  real  and  positive  limits  the  number  of  minima  to  one  and 
the  number  of  maxima  to  one  unless  V  4  ^  and  ;z|  satisfies  the  inequality 


1  4  v  1  -  4li/!® 
*  Z 


(244) 


When  v  is  complex  the  values  of  m  become  the  roots  of  a  quartic  equation. 

The  absolute  rounding  error  of  the  asymptotic  series  is  determined  by  the  relative 
rounding  error  in  the  largest  term  of  the  series,  while  the  absolute  truncation  error 
of  the  series  is  determined  by  the  magnitude  of  the  smallest  term  of  the  series  In 
order  to  control  the  relative  rounding  error  in  the  series  it  is  necessary  to  limit  the 
evaluation  of  the  series  to  conditions  where  the  terms  diminish  after  the  first  term 
The  full  accuracy  of  the  computer  is  guaranteed  only  if  the  order  and  the  argument 
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satisfy  the  criterion 


In  order  to  control  the  truncation  error  in  the  series  it  is  necessary  to  limit  the 
evaluation  of  the  series  to  conditions  where  the  smallest  term  is  on  the  order  of 
rounding  error.  For  the  CDC  6600  computer  the  argument  must  satisfy  the  criterion 
2  ;  <5  17.5. 

The  criteria  for  full  accuracy  severely  restrict  the  range  of  order  and  argument  in 
which  the  classical  series  may  be  applied. 

Required  for  the  Debye  approximation  are  polynomials  um(t )  which  are  expressed 
by  the  equation 

m 

Um(t)  =  E  cmt  (246) 

k  =  0 

for  which  the  coefficients  are  given  by  the  Amos6-10  recurrence  equations  The  evaluation 
of  the  coefficients  is  started  with  the  equation 

c00  -  1  (247) 


and  is  continued  with  the  equation 


°mk  =  I  ru  +  2k 


m  +  2k  1  f  | 

n  cm  - 1  ,t  "  ~  ~ 

c  [m  + 


m  -  Zk  -  3 


A  parameter  y  is  defined  by  the  equations 


tanh  v  =  v  1 

\  v 


1  -  >og - g  (250) 

1  +  V  ' 

Let  functions  S!  and  s2  be  defined  by  the  equations 

c+WUnh>-,)  »  U^COth?) 

S,  -  - p=5TT.----=-  r-rr  >  - - -  -  (cDl) 

\  yZ-nv  tanh  y  m-o  ( 1 1')  , 

e-rtUnh>.T,  X  u  (cothy)' 

s2  ,  • .  E  Am  (2  o2) 

\  2,-nv  tanh  y  m  o  (  L’l 

Thus  s2  is  obtained  from  s,  by  a  reversal  of  the  sign  of  n.  Both  functions  are  the 
product  of  a  factor  and  a  series.  The  factor  contains  two  radicals  and  a  logarithm, 
while  the  series  is  asymptotic.  Let  s,  be  the  Debye  approximation  with  a  positive  sign 
assigned  to  the  radical 


and  let  s2  be  the  Debye  approximation  with  a  negative  sign  assigned  to  this  radical 
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Reversal  of  the  sign  of  the  radical  replaces  the  argument  of  the  logarithm  by  its 
reciprocal  and  reverses  the  sign  of  the  logarithm.  Reversal  of  the  sign  of  the  radical 
replaces  the  exponential  function  by  its  reciprocal.  Thus  where  s,  is  large,  sa  is  small, 
and  vice  versa. 

The  radical  is  zero  at  branch  points  in  the  complex  n-plane  where  u  =  ±z.  In  the 
vicinity  of  each  branch  point  the  logarithm  may  be  expanded  in  powers  of  the  binomials 

1  -  2  or  1  +  -  (254) 

u  v 

Cancellation  of  the  lowest  order  terms  and  omission  of  the  highest  order  terms  lead 
to  the  approximation 


3 


(255) 


Three  nodal  lines  emanate  from  each  branch  point.  On  each  nodal  line, 

is,|  -  is2!  or  eta'ri/*s1j  -  s2  (256) 

For  the  exponential  functions  in  the  approximations  to  have  the  same  absolute  values 
it  is  necessary  for  their  arguments  to  be  pure  imaginary.  For  the  approximation  of 
the  argument  to  be  pure  imaginary  the  order  must  be  given  by  one  of  the  equations 


,  3  U 

u  z  -  A  u3  e3 

1/  i  Z  A  I/3  B3 


(v~  -  z)  (257) 
(v  ~  2)  (258) 


where  A  is  a  real  parameter.  Thus  the  nodal  lines  emanate  from  the  branch  points  in 
the  directions  of  the  three  roots  of  1  with  respect  to  a  line  which  makes  an  angle 
with  the  real  axis  equal  to  one  third  of  the  angle  which  z  makes  with  the  real  axis 
As  z  rotates  the  nodal  lines  rotate  at  one  third  the  rate  of  rotation  of  z. 

On  the  positive  side  of  the  imaginary  axis  computation  shows  th.it  the  Bessel  function 
is  given  by  S(  alone  along  the  nodal  lines  which  emanate  to  the  right  of  z,  but  the 
Bessel  function  is  given  by  s,  >  s2  along  the  nodal  line  which  emanates  to  the  left 
The  boundary  between  the  regions  where  s,  alone  is  used  and  where  s,  *  s2  must  be 
used  presumably  is  intermediate  between  nodal  lines,  where  s2  is  so  much  smaller 
than  s,  that  the  difference  in  formulation  is  immaterial. 

On  the  negative  side  of  the  imaginary  axis  computation  shows  that  the  Bessel 
function  is  given  by  the  sum 

s,  ‘  s2  4-  e42^,  (259) 

on  the  left  of  an  hyperbola  and  is  given  by  the  sum 

+  s,  (■  Sj  el2,,,'4s1  (260) 

on  the  right  of  the  hyperbola  In  these  formulae  the  ±  sign  is  -  in  a  region 
counterclockwise  from  an  outward  extension  of  z  and  is  +  in  a  region  clockwise  from 
the  extension  of  z  The  boundary  of  the  region  where  etZ"'l'ts1  must  be  included 
presumably  is  on  a  line  midway  between  nodal  lines,  where  etS’Tl''s1  is  so  much  smaller 
than  s,  4  s2  that  the  difference  in  formulation  is  immaterial. 
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For  a  large  value  of  v  on  the  positive  side  of  the  imaginary  axis  the  signs  of  radicals 
are  positive  and  the  logarithm  requires  no  correction,  but  in  other  regions  corrections 
are  required  because  the  radicals  and  logarithm  are  evaluated  on  the  assumption  that 
their  arguments  range  from  a  value  just  more  than  -rr  through  the  value  ~n.  The 
conditions  for  correction  are  determined  from  Argand  diagrams  in  the  complex  n-plane. 

Let  z  and  v  be  expressed  in  terms  of  their  components  by  the  equations 

2  =  r4-iy  v  =  A  +  i  y  (261) 


Let  A  diminish  from  +  <*>  with  the  other  components  constant.  Then  the  trace  of  v  is 
a  straight  line  parallel  to  the  real  axis  and  the  trace  of  1/n  is  a  circle  of  radius  1  2y 
which  starts  at  the  origin  and  has  a  center  offset  along  the  imaginary  axis  to  a  distance 
-1/2/u.  Multiplication  by  z  amplifies  and  rotates  the  circle.  The  ratio  z  v  is  a  circle 
which  intersects  the  origin  and  is  rotated  through  an  angle  equal  to  the  phase  of  z. 
The  trace  of  the  square  zz/us  is  a  cardioid  which  is  rotated  through  an  angle  equal 
to  twice  the  phase  of  z. 

The  sign  of  the  radical 


(262) 


is  positive  for  large  \i>\  but  must  be  reversed  whenever  its  argument  crosses  the 
negative  real  axis.  Let  the  argument  be  expressed  by  the  equation 

,2 

1  -5  rz  (263) 

v 

where  r  is  a  real  variable.  Solution  gives  the  equation 

-  i  \  1  -  r2  (26-1) 


For  the  ratio  z  V  to  be  real  it  is  necessary  for  the  components  to  be  related  in 
accordance  with  the  equation 


y  _  y 

x  A 

Then  the  ratio  z/v  is  given  by  the  equation 


z 

V 


y 

M 


\  1 


(265) 


(266) 


The  argument  crosses  the  negative  real  axis  therefore  on  the  line  where  v  is  congruent 
to  z  and  y  y 

The  argument  of  the  logarithm  is 


2 

V 


(267) 


Even  when  'in.  is  large,  the  argument  crosses  the  negative  real  axis  when  y  and  y  arc 
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opposite  in  sign.  Let  the  argument  be  expressed  by  the  equation 

Z 

V  1 


(268) 


where  r  is  a  positive  real  value.  Solution  gives  the  equation 

2  _  2r 
v  1  +  r2 


(269) 


This  function  of  r  ranges  from  zero  through  a  minimum  and  back  to  zero  as  r  ranges 
from  0  to  “.  The  minimum  is  located  where  the  derivative  with  respect  to  r  is  zero. 
The  minimum  value  is  - 1 .  For  the  ratio  z/v  to  be  real  it  is  necessary  for  the  components 
to  be  related  in  accordance  with  the  equation 


y  =  M 

x  \ 

Then  the  ratio  z/v  is  given  by  the  equation 

2  y  2  r 

v  m  1  +  r2 


(270) 


(271) 


The  argument  crosses  the  negative  real  axis  therefore  on  the  line  where  u  is  congruent 
to  2  and  /z,  y  have  opposite  signs  with  \fi\  >  jy|.  At  the  crossing  the  logarithm  is  corrected 
by  *2ni. 

The  sign  of  the  radical 


1 

{vz~  z2)* 


(272) 


is  positive  for  large  \u\  but  must  be  reversed  whenever  its  argument  crosses  the 
negative  real  axis.  Let  the  argument  Le  expressed  by  the  equation 

i/2  -  22  =  -  r2  (273) 

where  r  is  a  real  variable.  For  the  argument  to  be  real  it  is  necessary  for  the  components 
to  be  related  in  accordance  with  the  equation 

V  =  xy  (274) 


Then  the  argument 


z 2  is  given  by  the  equation 


2  z  *zyz 

vi  -  z2  -  v- 

M 


x2  +  y 2 


(275) 


The  argument  is  positive  only  if  |/u|  <  |y|  whence  its  square  root  can  be  real  The  trace 
of  the  condition  for  correction  is  an  hyperbola  with  the  equation 


\  = 


xy 


(276) 
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Whether  a  correction  is  necessary  depends  upon  the  sign  of  the  radical 


The  correction  is  necessary  only  if  the  radical  has  a  direction  opposite  to  the  complex 
conjugate  of  v.  The  condition  for  the  correction  is  not  met  when  the  sign  of  the  radical 
has  been  reversed. 

The  Debye  approximation  is  not  useful  inside  a  boundary  where  the  smallest  term 
of  the  series  would  be  greater  than  rounding  error.  Exploratory  computations  have 
shown  that  the  limiting  boundary  for  large  argument  and  order  is  a  cubic  parabola. 
An  empirical  criterion  which  defines  the  limiting  boundary  for  accurate  evaluation  is 
given  by  the  equation 

!  1/  3  u  I3 

|  1  -  -  1  +  -  j 

,4|z|e— El  (278) 

where  the  constant  A  ~  0.004.  This  criterion  matches  the  true  boundary  for  extreme 
arguments  and  does  not  deviate  significantly  from  the  true  boundary  at  intermediate 
arguments.  If  the  criterion  is  not  met  for  a  given  value  of  v,  then  unit  increments 
are  added  to  u  until  the  criterion  is  met. 

Inasmuch  as  each  term  of  the  Debye  series  is  itself  a  polynomial  it  has  a  number 
of  roots  equal  to  its  degree.  In  order  to  avoid  a  premature  termination  of  the  series, 
the  polynomials  are  evaluated  both  with  the  complex  arguments  and  with  their  absolute 
values  until  the  summation  of  the  absolute  values  of  the  terms  remains  stationary. 

Values  of  the  order  and  the  argument  which  are  outside  the  range  of  the  classical 
series  still  can  be  reached  by  an  application  of  recurrence  equations  which  are  started 
within  the  range  of  the  classical  series.  The  Bessel  functions  and  the  Weber  functions 
satisfy  the  recurrence  equations 

Zv 

Vi(*)+^«(*)=  —  -M*)  (279) 


i(z)  +  *V+i(z)  =  —  >V(z) 


and  the  recurrence  relation 


■Mz)}'„n(z)  ~  ^k+i(z)>'v(z)  = -  (281) 

7TZ 

Let  eu  be  the  rounding  error  which  has  been  introduced  in  the  ^tth  cycle  of  iteration 
The  persisting  error  in  the  yth  cycle  is  given  by  the  expression 


| Y^t(z)Jv(z)  -  ^*,(z)}'i/(z)j 


for  a  descending  recurrence,  and  by  the  expression 


+  g z  jlViUVxU)  -  i(z)ru(2)j 
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for  an  ascending  recurrence  In  the  computation  of  Jv(z)  the  recurrence  must  be 
cycled  in  whichever  direction  !)r„(z)|  diminishes  relative  to  ,Jv(z)\. 

The  recurrence  is  applied  actually  to  r(v  +  \){^z)~vJv(z)  when  the  convergent  series 
is  the  origin  of  Jv(z)  in  order  to  keep  the  recurrence  within  the  index  range  of  the 
computer.  The  recurrence  is  applied  actually  to  (ttz/2 )i/aJ„(z)  when  the  asymptotic 
series  is  the  origin  of  J„(z)  in  order  to  improve  the  efficiency  of  computation. 

The  convergent  approximation  is  used  when  z\  g  17.5.  There  are  two  modes  of 
computation  with  the  convergent  approximation.  In  the  first  mode  the  convergent 
series  is  used  for  orders  with  large  positive  real  parts,  and  descending  recurrence  is 
used  to  bring  the  order  down  to  orders  with  less  positive  real  parts.  In  the  second 
mode  the  convergent  scries  is  used  for  orders  with  large  negative  real  parts,  and 
ascending  recurrence  is  used  to  bring  the  order  up  to  orders  with  less  negative  real 
parts.  The  relative  error  in  either  recurrence  increases  until  the  order  crosses  a  nodal 
line,  then  the  error  remains  constant  thereafter. 

Boundaries  between  modes  are  located  where  the  errors  are  the  same  for  both 
modes.  The  ideal  boundaries  between  modes  are  represented  by  complicated  surfaces 
in  the  four-dimensional  space  of  order  and  argument.  Information  about  the  location 
of  boundaries  is  derived  from  comparisons  between  single-precision  and 
double-precision  computations.  The  boundaries  can  be  perceived  only  dimly  in  the 
computations  because  of  random  fluctuations  in  rounding  error  Within  the  random 
fluctuations  the  boundaries  can  be  simulated  by  surfaces  with  polygonal  sections. 

In  the  first  mode  of  computation,  recurrence  is  started  with  that  order  with  a 
positive  real  part  which  satisfies  the  criterion 

\v\  >  Ijzi*  (284) 

In  the  second  mode  of  computation,  recurrence  is  started  with  that  order  with  a 
negative  real  part  which  satisfies  the  criterion 

Xe  v  b  KM  -  '3m  i> )  (285) 


The  first  mode  is  used  in  preference  to  the  second  mode  when  the  order  and  argument 
satisfy  both  of  the  criteria 


and 

Re  v  -  ’(  z 

-  \.3m  v 

3m  z  ^  3m  v  ) 

(286) 

3te  v  4  z  B  5  z 

3m  v  ■  3m 

z  as'”  3m  v  -  3m  z‘ 

(287) 

The  two  criteria 
accordance  with 

combine  to  give 
the  computations. 

polygonal 

sections  on  planes  of  constant 

z  in 

The  Debye  approximation  is  used  when  z  ■  17  5.  A  zone  in  the  v  plane  from  which 
the  Debye  approximation  is  excluded  can  be  traversed  by  descending  recurrence  when 
the  zone  is  on  the  positive  side  of  the  imaginary  axis  Only  part  of  the  zone  can  be 
traversed  by  descending  recurrence  when  the  zone  is  on  the  negative  side  of  the 
imaginary  axis  An  ascending  recurrence  with  the  Debye  approximation  gives  too  much 
error  Within  the  /.one  of  exclusion  it  is  possible  to  use  the  convergent  approximation 
Insofar  as  the  descending  recurrence  starts  with  equally  accurate  initial  values  for 
either  the  Debye  approximation  or  the  convergent  approximation,  the  error  during 
descending  recurrence  is  Ihe  same  for  both  approximations  The  same  boundaries 
apply  between  descending  recurrence  from  the  Debye  approximation  and  ascending 
recurrence  from  the  convergent  approximation  The  Bessel  function  is  computed  with 
a  combination  of  Debye  approximation  with  descending  recurrence  and  convergent 
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approximation  with  ascending  recurrence 

The  relative  error  is  not  uniform  over  the  complex  u-plane.  Wherever  the  Bessel 
function  approaches  zero  the  relative  error  approaches  infinity.  The  relative  error  is 
large  over  a  nodal  line  where  the  Bessel  function  is  small.  When  z  is  real,  the  Bessel 
function  can  be  computed  with  full  machine  accuracy.  When  z  is  rotated  out  of  the 
real  axis  a  zone  of  rounding  error  appears  and  grows.  The  zone  of  rounding  error 
straddles  the  negative  real  axis.  Eventually  the  zone  of  exclusion  cuts  off  the  zone  of 
rounding  error. 

Programming 

SUBROUTINE  CBSSlJ  (AZ.  CN,  Fj) 

*%#**#*******#*****%*******.%****%*****.  t**4*«*«**«x**«******#***«***«**«***»:**«*r*  ****4*****»»* 

FORTRAN  SUBROUTINE  FOR  COMPLEX  BESSEL  FUNCTION  OF  FIRST  KIND 

*n*iH[***m.*#  +  *#**#mc  +  *********m  +  **iti*m*#*****  +  **********m*****  +  ******m*mm***4:»******»**» 

The  real  and  imaginary  parts  of  the  complex  argument  z  are  given  in  array  ~Z, 
and  the  real  and  imaginary  parts  of  the  complex  order  u  are  given  in  array  C\.  The 
complex  Bessel  function  is  computed  with  series  expansions  and  recurrence  relations 
Calls  are  made  to  Subroutine  CGAMMA.  The  real  and  imaginary  parts  of  the  complex 
function  Jv{z)  are  stored  in  array  rJ. 

SUBROUTINE  DBSSlJ  (AZ,  CN,  FJ) 

******#************#**************«**********«*•****•««#*****#*+»**********««****•*••»•****»• 

FORTRAN  SUBROUTINE  FOR  DOUBLE  -  PRECISION  BESSEL  FUNCTION  OF  FIRST  KIND 
****«****«**«**«*****««*«***«********«**.********«•**-******«***•***•*«*«************•*«•*»•«•* 

The  real  and  imaginary  parts  of  the  double  precision  argument  z  are  given  in  array 
AZ,  and  the  real  and  imaginary  parts  of  the  double  precision  order  v  are  given  in 
array  CL.  The  double-precision  Bessel  function  is  computed  with  series  expansions 
and  recurrence  relations.  Calls  are  made  to  Subroutine  DG-'.'VA.  The  real  and  imaginary 
parts  of  the  double  precision  function  J„(z)  are  stored  in  array 


DISCUSSION 

The  incomplete  beta  function  B(p,  q,  x)  is  defined  by  the  equation 

B(p,  q.x)  f  fp  ‘(1  -  <)«'  1  dt  (288) 

Jo 

and  the  incomplete  gamma  function  \'(p,  x)  is  defined  by  the  equation 

T(p,  x)  f  fp  1  e  ‘  df  (289) 

Jo 

A  new  subroutine  BL'AX  computes  directly  the  incomplete  beta  function  and  a  new 
subroutine  GAMMAx  computes  directly  the  incomplete  gamma  function  for  arbitrary 
real  order  and  real  argument  More  efficient  subroutines  for  the  incomplete  beta  ratio 
of  half  integer  order  and  the  incomplete  gamma  ratio  of  arbitrary  order  have  been 
prepared  by  A  R  DiDonato25,28 

The  potential  of  a  rectangular  plate  is  useful  for  the  computation  of  flow  around 
struts27  and  hulls28 

The  single  precision  Bessel  functions  wore  checked  by  comparison  with 
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double  precision  computations  The  single  precision  Jn{z)  was  compared  directly  with 
the  double  precision  Jn(z).  The  single- precision  Y „(z)  was  compared  with 
double  precision  values  from  the  equation 


irau 

Tt  di/ 


where  the  partial  derivatives  with  respect  to  v  were  estimated  by  finite  differences 
between  values  with  v  n  t  c.  The  single  precision  In(z)  was  compared  with 
double  precision  values  from  the  equation 

/„(*)  =  (  -i)Vn(i5)  (291) 

when  y  was  negative  and  from  the  equation 

In(z)  ---  inJn(  iz)  (292) 

when  y  was  positive.  The  single  precision  A n(z)  was  compared  with  double  precision 
values  from  the  equation 


L  dl'  Un 


where  the  partial  derivatives  with  respect  to  u  were  estimated  by  finite  differences 
between  values  with  v  n  t  e.  The  values  of  the  functions  Iv(z)  were  derived  from  the 
equation 

Uz)  -  e  v*iJl.(iz)  (29-1) 

when  y  was  negative,  and  from  the  equation 

l As)  ev*%(--iz)  (295) 

when  y  was  positive.  The  first  order  difference  across  a  point  is  in  error  only  in  the 
third  order,  and  the  accuracy  of  the  difference  is  one  and  a  half  precision  when  the 
difference  in  argument  is  at  the  half  precision  level. 

On  page  265  of  Theory  of  Bessel  Functions  by  Watson2  there  is  a  figure  which  gives 
the  boundaries  for  various  combinations  of  st  and  s2  in  the  u  z-plane.  This  figure 
agrees  with  the  analysis  herewith  for  real  z  but  does  not  show  the  i  station  of  boundaries 
at  one  third  the  angle  of  rotation  of  z  which  is  characteristic  of  the  nodal  lines 
Subroutines  for  Jv(z )  and  Iw(z)  have  been  programmed  by  Amos.  Daniel,  and  Weston10 
Their  subroutines  are  valid  for  positive  real  argument  and  positive  real  order  They 
supplement  the  classical  series  with  the  Debye  approximation,  and  they  use  the  Olver 
approximation  where  the  Debye  approximation  is  ineffective  A  subroutine  for  1  \(z) 
has  been  programmed  by  Cody,  Motley,  and  F’ullerton u.  Their  subroutine  is  valid  for 
positive  real  argument  and  positive  real  order.  It  uses  a  sequence  of  Taylor  scries 
expansions  for  lowest  orders,  and  it  uses  ascending  recurrence  to  reach  higher  orders 
The  names  of  the  subroutines  for  all  four  Bessel  functions  are 


The  first  three  subroutines  have  been  checked  against  double  precision  computations 
Four  subroutines  in  the  present  project  are  valid  for  complex  argument  but  for 
integer  order  They  supplement  the  classical  series  with  rational  approximations  at 
small  orders  and  use  recurrence  to  extend  the  range  of  orders  The  names  of  these 


subroutines  are 


BSSLJ  BSSLY  BSSLI  BSSLK 

The  new  subroutines  are  more  compact  than  the  other  subroutines.  They  are  more 
efficient  where  they  use  rational  approximations.  Otherwise  the  new  subroutines  are 
as  efficient  at  low  orders  as  the  other  subroutines  would  be  if  they  were  converted 
to  complex  arithmetic.  Conversion  of  arithmetic  from  real  to  complex  may  increase 
the  time  of  computation  by  1§  on  the  CDC  6600  where  there  can  be  parallel  processing. 
The  accuracies  of  the  new  subroutines  for  low  orders  are  within  one  digit  of  the 
accuracies  of  the  other  subroutines.  The  efficiencies  and  the  accuracies  of  the  new 
subroutines  deteriorate  with  increase  of  order  into  the  range  of  order  where  the  Debye 
approximation  is  used. 


CONCLUSION 

The  Debye  approximation  can  be  used  for  the  computation  of  Bessel  functions  of 
complex  order  and  complex  argument  everywhere  except  in  a  zone  of  rounding  error 
where  the  value  of  the  Bessel  function  is  relatively  small  In  the  zone  of  rounding 
error  the  convergent  series  with  ascending  recurrence  gives  better  accuracy  than  the 
Debye  approximation  with  descending  recurrence 
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APPENDIX  A 


LEGENDRE  FUNCTIONS 


LEGENDRE  FUNCTIONS 


An  analysis  of  Legendre  functions  of  integral  order  is  to  be  found  in  Modem  Analysis 
by  Whittaker  and  Watson1.  The  differential  equation  for  the  Legendre  functions  is 


(1 


_  dzw  dm  , 

zz) - ~  -  2 2  —  +  n(n  +  l)ic  =  0 

dz  dz 


The  Legendre  function  of  the  first  kind  is  given  by  the  Sehlafli  integral, 

1  f  (t2  l)n 


(i) 


(2) 


where  the  contour  of  integration  encircles  2  once  counterclockwise.  That  this  integral 
satisfies  the  differential  equation  is  easily  verified  since  substitution  in  the  differential 
equation  leads  to  the  circuit  integral, 


,  .  dzPn  dPn  .  n  ■*  1  _  ,  _ 

(1-s8)-- 2z  -  —  ”  +  n(n  -  1  )Pn  -  •  (p  ■■ -  3  a -  ~  .  >  dt 

v  '  dz*  dz  2  ‘m  J  dt  (  (t  z)  2  ) 


d  {(**-  l)n*‘) 


(3) 


That  the  integral  coincides  with  the  Legendre  polynomial  is  verified  when  the  integral 
is  evaluated  for  2  --  ±  1. 

The  Sehlafli  integral  may  be  derived  by  n  fold  differentiation  of  the  Cauchy  integral 


1  if  'ft, 

2rri  J  (t  z) 

Thus  the  Legendre  polynomials  are  given  by  the  Rodrigues  formula, 


P  (2)  . 1  rfn  (j2  l)n 

2nn!  dzn  K 


Expansion  of  the  identity 


1  f  d  J  (t*  1)"*'  ) 


2n‘am  J  dt  (  (t  2 )n * 1  ) 


dt  0 


leads  to  the  equation 

PnM(2)  2  PJZ) 
Expansion  of  the  identity 


I 


f  u2  nn 


dt 


\  f  l)nldt  0 

2 n-'m  J  dt  (  " 


leads  to  the  equation 


nzPJz)  n  Pn  ,U) 


(t  2)"  ) 


n  *  1 
2n',7TZ  J 


u2  ir 


dt 


(t  zV 

Elimination  leads  to  the  recurrence  equation 

nPn  ,(2)  (2n  •  1)2 /’„(  =  )  •  (n  •  \)Pn.,(z)  0 


’I) 


(5) 


(P) 


(T) 


(Hi 


(9) 


1 10) 
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and  differentiation  leads  to  the  recurrence  equation 

(2n+ i)pn(*)  =  p;M(z) (11) 

Thus  all  of  the  Legendre  functions  and  their  derivatives  may  be  synthesized  from  the 
two  functions  of  lowest  order. 

Differentiation  of  the  recurrence  equations  leaves  one  term  of  lower  order  in  each. 
Repeated  differentiation  with  elimination  of  the  term  of  lower  order  leads  to  the 
recurrence  equation 


(n  -  m) 


dmPn 


(2  n  +  1): 


dmP  dmP 

T  +  (n  -m  +  1 )  —-—Odl 

dz  '  '  dzm 


(12) 


which  connects  the  derivatives  of  functions  of  progressively  increasing  order.  Repeated 
differentiation  uf  the  differential  equation  leads  to  the  recurrence  equation 


(n  -  m)(n 


m  + 


1) 


dm  lPn 
dzm  1 


-  2mz 


d^P_n 

dzm 


(1  -z2) 


dm*iPn 

dz™'C 


=  0 


(13) 


which  connects  the  progressively  higher  derivatives  of  functions  of  common  order. 
The  Legendre  function  of  the  second  kind  is  given  by  the  integral. 


«nU) 


i  rl  (i  -t*)n 
2"*1  J  .  ,  (z  -  f)""1 


(14) 


Substitution  of  the  integral  in  the  differential  equation  and  in  the  recurrence  equations 
shows  that  Qn(z)  satisfies  the  same  equations  as  Pn(z)  provided  that  9le(n  ■*  1)  -0. 


2 


APPENDIX  B 


MULTIPLE-PRECISION  COMPUTATIONS 


MULTIPLE- PRECISION  COMPUTATIONS 
Formula  Algebraic  Processor 

FLAP,  or  Formula  Algebraic  Processor,  is  a  computational  system  which  has  been 
developed  by  Morris18  for  the  manipulation  of  mathematical  expressions.  The  FLAP 
system  refers  to  LISP,  or  List  Processor,  which  manipulates  lists  of  numbers  or  symbols. 
As  a  checkout  of  the  numerical  feature  of  FLAP,  high- precision  values  of  Bernoulli 
numbers  and  Riemann  zeta  functions  have  been  computed  Values  were  computed  to 
80  digits  with  an  accuracy  of  75  digits,  and  then  the  values  have  been  truncated  back 
to  64  digits.  The  values  have  been  punched  on  cards  and  are  printed  in  Appendix  C. 

Bernoulli  Numbers 

There  are  two  principal  versions  of  the  Bernoulli  numbers.  The  older  version  is 
derived  from  the  definition 


|zcot|z  =  1  -  £  Bn 

U=  1 


(2  n)\ 


(1) 


where  Bn  is  the  nth  Bernoulli  number  in  the  older  definition.  The  newer  version  is 
derived  from  the  definition 


—  'P  R  - 

-  1  "  n! 


(2) 


where  Bn  is  the  nth  Bernoulli  number  in  the  newer  definition.  The  Bernoulli  number 
Bn  in  the  older  definition  is  equivalent  to  {  -  \)n~  lBZn  in  the  newer  definition  The 
change  of  definition  adds  one  number  to  the  set  and  makes  some  formulae  more 
compact. 

Values  of  Bn  in  the  newer  definition  can  be  obtained  by  the  division  of  z  by  the 
series  expansion  for  e*-  1.  The  first  few  values  are  given  by  the  equations 

Ba  =  1.  B,  =  -J.  Bz  =  l  B3  -  0,  B<  ----  -±0,  0,  Be  ^  (3) 


If  z  is  replaced  by  -z  in  the  definition  equation  to  give  terms  which  are  subtracted 
from  the  terms  of  the  definition  equation,  then  cancellation  of  terms  leads  to  the 
equation 


2  V'  /?„  _ 

■  2n  *  i  /o  1  \  l 

n.,o  (2n  ‘  ')! 


(4) 


Comparison  of  coefficients  of  powers  of  z  shows  that  the  Bernoulli  numbers  of  odd 
order  satisfy  the  equation 


Bint  i  -  0 


(n  ■  0)  (5) 


If  z  is  replaced  by  ~z  in  the  definition  equation  to  give  terms  which  are  added  to  the 
terms  of  the  definition  equation,  then  cancellation  of  terms  leads  to  the  equation 


\z  coth  |z 


K  BtU  (2 nj! 


(6) 
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Substitution  of  2 iz  for  z  leads  to  the  equation 

d,.  i 

„o  (2~)' 

which  gives  the  Maclaurin  expansion  for  the  cotangent.  From  the  identity 

2  cot  2z  =  cot  z  -  tan  z 


(V) 

(6) 


may  be  derived  the  equation 


tan  z  = 


(-  l)n22n(22n 


E 

n=l 


1  )BZ 


(2  n)! 


(9) 


which  gives  the  Maclaurin  expansion  for  the  tangent. 

If  z  is  real,  then  integration  around  a  contour  which  is  bounded  by  the  real  axis, 
a  line  at  i  parallel  to  the  real  axis,  and  the  imaginary  axis,  but  does  not  contain  the 
points  0,  i,  leads  to  the  equation 

sin  zt  .  .  1  ,  , 

,2wt  _~7  dt  =  4  coth  Yz  (i.- real)  (10) 

Substitution  for  coth|z  its  series  expansion  in  terms  of  Bernoulli  numbers,  substitution 
for  sinzf  its  series  expansion  in  powers  of  zt,  and  comparison  of  coefficients  shows 
that  the  even-ordered  Bernoulli  numbers  are  given  by  the  equation 

^2n  1 

fl2n  =  -(-Dn4n  -z^rr-dt  (11) 

J  o  e  1 

A  change  of  variable  in  the  integrand  expresses  the  Bernoulli  numbers  in  terms  of 
Riemann  zeta  functions  of  even  order. 

Values  of  Bernoulli  numbers  have  been  obtained  from  Tables  of  the  Higher 
Mathematical  Functions  by  Davis23.  They  are  listed  in  Table  1  at  the  same  level  of 
accuracy  as  the  other  numbers  which  have  been  obtained  through  FLAP.  Values  of 
Bernoulli  numbers  for  use  in  the  computer  were  computed  by  the  method  of  Knuth 
and  Buckholtz22.  The  derivative  of  the  nth  power  of  the  tangent  is  given  by  the  equation 


—  tannz  =  n  tan"  'z(l  +  tan2z) 
dz 


(12) 


Thus  the  derivative  of  a  polynomial  in  powers  of  tan  z  remains  a  polynomial  in  powers 
of  tan  z.  The  derivative  of  the  power  polynomial 

,v 

E  an  tannz  ( 1 3) 

n-0 

is  obtained  by  the  transformation 

an  -*(**-  l)a„_,  +  (n  +  l)an»,  (14) 

If  a,  is  unity  and  all  other  coefficients  an  are  zero  initially,  then  iteration  of  the 
transformation  gives  all  of  the  derivatives  of  tan  z.  from  which  numerical  values  for 
the  coefficients  of  the  Maclaurin  series  can  be  computed  From  the  coefficients  of  the 
Maclaurin  series  are  computed  the  Bernoulli  numbers 


1 


2 


Bernoulli  Polynomials 

Bernoulli  polynomials  were  derived  from  an  unconventional  generating  function  by 
Whittaker  and  Watson1  Their  polynomials  were  terminated  with  terms  in  z  or  c2, 
whereas  the  conventional  polynomials  are  terminated  with  terms  in  1  or  2.  The 
conventional  definition  is  to  be  preferred  because  then  the  polynomials  are  terminated 
automatically  when  factorials  of  negative  integers  appear  in  the  denominators  of  their 
terms. 

Bernoulli  polynomials  are  defined  by  the  equation 


f  eTTT  =  E  B-(z}  n7  (15> 

1  n=0  n- 

where  Bn( z)  is  the  nth  Bernoulli  polynomial.  Term  by  term  multiplication  of  the  series 
expansion  for  t/(e‘  -  l)  by  the  scries  expansion  for  ezt,  and  comparison  of  coefficients 
leads  to  the  equation 

»■*"'*  (,6» 

Thus  the  nth  polynomial  Bn(z)  is  of  the  nth  degree  Differentiation  of  the  terms  of 

this  polynomial  n  -  m  times  is  expressed  by  the  equation 

m.  | 

5^®  1,71 

Special  values  for  2  =  0  are  given  by  the  equations 

Bn(0)  =  Bn  (m  n)  ( 18) 


Sin-m)(0)  ---  B„ 

m  I 


(m  Ml  (10) 


Thus  the  Bernoulli  numbers  are  the  constants  in  the  Bernoulli  polynomials 

Subtraction  of  the  terms  in  the  generating  equation  from  the  same  terms  with  2 
replaced  by  z  +  1  leads  to  the  equation 

i  \Bn(z+  1)  Bn(z)\  1  (d0) 

n=0  1  )  n' 

Series  expansion  of  e11  and  comparison  of  coefficients  lead  to  the  equation 

Sn(2  +  1)  -  Bn(z)  -  n:n  1  (SI) 

Differentiation  of  the  terms  of  this  equation  n  m  times  is  expressed  by  the  equation 

a<r  m)(2  +1)  B{nn  m)(z)  -  -  n!  -  2m  1  (22) 

(m-  1)! 

Special  values  for  z  0  are  given  by  the  equations 

ain  l)(l)  Bin  u(0)  n!  (m  1)  (S3) 


a<T  m)(l)  B{”  m'(0)  0 


(m  1)  (21) 


3 


Thus  the  Bernoulli  polynomials  are  the  same  at  both  ends  of  the  range  0  ^  2  S  1  if  the 
order  n  *  1. 

A  comparison  of  coefficients  of  the  powers  of  t  in  the  identity 

e-*t 


shows  that  the  Bernoulli  polynomials  satisfy  the  equation 

Bn(l  -z)  -  (  l)n£nU) 


(25) 

(26) 


Thus  the  even-ordered  polynomials  are  symmetric  and  the  odd-ordered  polynomials 
are  antisymmetric  with  respect  to  the  midpoint  of  the  range  0  S  2  S  1. 

A  comparison  of  coefficients  of  the  powers  of  t  in  the  identity 


t 


(27) 


shows  that  special  values  for  2  -  |  are  given  by  the  equation 

Bn{\)  -  (2l  n  •  1  )Bn  (28) 

Thus  the  Bernoulli  polynomials  have  opposite  signs  at  the  midpoint  and  at  the  ends 
of  the  range  0  5?  2  §  1. 

Differentiation  of  the  Bernoulli  polynomials  is  expressed  by  the  equation 

B'n(z)  nBn.i(z)  (29) 

Integration  of  the  Bernoulli  polynomials  is  expressed  by  the  equation 


Thus  the  integral  of  a  Bernoulli  polynomial  is  zero  in  the  range  0  l-  z  r-  1  if  the  order 
n  *  0. 

Further  differentiation  of  the  Bernoulli  polynomials  is  expressed  by  the  equation 

B'n(z)  -  n(n  l)Bn.a(z)  (31) 

The  polynomial  B,(z)  of  first  order  is  monotonic  with  a  uniform  sign  throughout  each 
side  of  the  midpoint  of  the  range  0  s?  2  ^  1.  The  second  derivative  of  the  next  higher 
odd  ordered  polynomial  has  a  uniform  sign  throughout  either  half  of  the  range 
0  S  2  g  1,  and  the  polynomial  has  roots  only  at  0.  1  By  induction  it  may  be  inferred 

that  all  higher  odd-ordered  polynomials  have  the  same  property  The  first  derivative 
of  each  even-ordered  polynomial  has  a  uniform  sign  throughout  each  side  of  the 
midpoint  of  the  range  0  %  z  ^  1.  Each  even  ordered  polynomial  has  only  two  roots  in 
the  range  0  z  S;  1.  The  absolute  value  of  every  even  ordered  polynomial  is  a  maximum 
at  the  ends  of  the  range 

Euler-Maclaurin  Expansion 

An  elegant  derivation  of  the  Euler  Maclaurin  formula  has  been  given  by  Whittaker 
and  Watson1.  However,  their  unconventional  definition  of  the  Bernoulli  polynomials 
warrants  a  review  of  their  derivation  They  started  with  a  formula  which  they  attributed 
to  Darboux. 

Let  f(z)  be  analytic  at  all  points  on  a  straight  line  which  joins  a  to  2,  and  let  ^(t) 


4 


OirWM*  - 


be  any  polynomial  of  degree  n  in  t.  Then  if  0  t  §  1.  differentiation  is  expressed  by 
the  equation 

H  u 

I'  (  - 1 ) m( z  -  a)"V<n~m)(£)/<m,(a  +  (z-  a)t ) 

dt  m-  1 
n 

■=  E  (  -l)m(2  -  a)m4«<n-m"I)(f)/<m)(a  +  (*  -  a)f) 

m  -  l 

+  E  (  l)m(z  -  a)m+Vn~m)(f)/<m+1>(a  +  (2  -  a)t)  (32) 

m=  l 

In  the  first  summation  on  the  right  side  of  the  equation  the  substitution  of  m  +  1  for 
m  leads  to  the  cancellation  of  all  terms  of  the  summations  except  the  first  term  of 
the  first  summation  and  the  last  term  of  the  second  summation.  Thus  the  equation 
can  be  reduced  to 

-%  E  (-l)m(2  -  a)V,l'm)(f)/(m)(a  +  (z  -  a)t) 

dt  m=  1 

=  -  (z  -  a)<fi{n\t)f'(a  +  (z  -  a)f)  +  (-  1)"(2  -  a)n*l<fi(t)f(n'  l,(a  +  (z  -  a)<)  (33) 

Inasmuch  as  is  a  polynomial  of  nth  degree,  the  nth  derivative  is  constant 

and  equal  to  y<nl(0)  or  ^<n)0)-  Integration  with  respect  to  t  from  0  to  1  gives  the 
equation 


V  (-  l)m(2  -  a)m 

m=0 

=  (-  l)n(2  -  a)ntl  f‘  <fi(t)fn*l)(a  +  (2  -  a)f)  dt  (34) 

Jo 

Let  the  polynomial  ip(t)  be  identified  with  the  polynomial  Bn(t). 

Substitution  of  the  special  values  for  the  two  Bernoulli  numbers  of  lowest  order  in 
the  Darboux  identity,  and  substitutions  of  2 k  for  m  and  2n  for  n  lead  to  the 
Euler-Maclaurin  expansion13  15 

f{z)  -  f(a)  =  -g(z  -  a)[f(z)  +  /'(“)) 


*><n~m,(l)/(m)(2)  -  V.(n-m>(0)/«m>(a) 


f  (z  ~  a)zkBzk 
(2fc)! 


| f(ik)(z)  ~f{zk)( a) 


+ 


(z  -  a)2n+1  cl 

(2n),  -  J  San(f)/(en+,)(a  +  (2  -  a)f)  dt 


(35) 


Integration  by  parts  in  this  equation  continues  the  summation  through  additional 
terms. 

Addition  of  the  nth  term  of  the  summation  to  the  integral  in  the  equation  is 
equivalent  to  the  replacement  of  BZn(t)  in  the  integrand  by  BZn(t)  -  BZn.  Let  the 
difference  Bgn(t)  -  Bgn  be  replaced  by  its  average  value  -Bgn.  Then  the  expansion  is 
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expressed  by  the  equation 

/(*)  -  /(a)  =  !(z  -  a)(/'(z)  +/'(a)t 


y  (2  - 

...»  (2fc)! 


!/<2*>(2)-/<2t>(a) 


(2  -  a)in*1B 


(2n)i 


—  /<2n+1)(a  +  (2  -  a)0) 


(36) 


where  0  is  some  number  in  the  range  0S6J  1.  The  first  term  on  the  right  side  of 
this  equation  expresses  the  trapezoidal  rule  for  /  f'(t)  dt,  the  summation  provides  a 
correction  for  the  trapezoidal  rule,  and  the  last“term  provides  an  estimate  of  the 
error  in  the  correction. 

Riemann  Zeta  Function 

The  Riemann  zeta  function  {(s)  is  defined  in  terms  of  its  argument  s  by  the  equations 


f(s) 


“  1  _  1  p  x»-* 

Si  **  r(s)  J0  e1  -  1 


(37) 


Direct  evaluation  of  the  series  to  an  acceptable  level  of  accuracy  is  not  feasible  for 
a  small  value  of  s.  In  the  special  case  where  s  is  an  even  integer  2n,  the  Riemann 
zeta  function  is  given  by  the  equation 

(-l)B(2n  )ZnBZn 


~  2(2n)! 

Otherwise  application  of  the  Euler-Maclaurin  expansion  gives  the  equation 


(36) 


<T(s)  =  E  r. 


1 


1 


T,  k*  2K’  (s  DA"*'1 


+  E 


r(s) 


1 


l=l  (2m)!  r(s  +  2m  -  2)  A',*2m  1 
Bzm  r(s)  "  1 


V 


(2 M)\  r(s  +  2M  -  1)  £k  {k  +  0)’* 


(39) 


Consideration  of  the  magnitude  of  the  remainder  by  Morris  has  indicated  that  if  A'  41 
and  si’ 2,  then  M  ^  39  for  an  error  of  less  than  4xl0~74.  High-precision  values  of  the 
Riemann  zeta  function  for  integer  values  of  the  argument  are  listed  in  Table  II  The 
most  accurate  values  which  previously  have  been  computed  only  had  fifty  digits21 

Gamma  Function 

The  logarithm  of  the  gamma  function  is  defined  by  the  equation 


log  T(  1  +  2)  =  -  yz  +  £  ]  -  +  log  n  -  log(n  +  z)  j 


(40) 


The  first  derivative  is  given  by  the  equation 

d 


dz 


log  r(t  +  z) 


y  - !_j 

1  (  n  n  +  z 


6 


(41) 


(42) 


i: 


If 


‘i 


The  higher  derivatives  are  given  by  the  equation 

57* log  r(1  +  2)  =  (~1)m(m ' 1)!  £  (^r 

Thus  the  Taylor  series  expansion  of  the  logarithm  of  the  gamma  function  is  given  by 
the  equation 

log  r(l  T  z)  =  -  yz  +  E  m)  zm  (43) 

m=2  m 

Addition  of  the  Taylor  series  expansion  for  log(l  +  z)  gives  the  equation 

log  T(2  +  z)  =  (l  -  y)z  +  E  -Um  (44) 

m=2  m 

Let  log  r(2  +  z)  be  expressed  by  a  series  in  accordance  with  the  equation 

00 

log  r(2  +  z)  =  -  E  atz*  (45) 

fc=l 

and  let  1/T(2  +  z)  be  expressed  by  a  series  in  accordance  with  the  equation 

~~  =  E  bmzm  (46) 

‘(2  +  z)  m,0 

Differentiation  and  substitution  gives  the  equation 

oo  oo  oo 

E  (n  *  IKm2"  =  E  bmZm  E  (fc  +  l)ai+iz‘  (47) 

n=G  m  =  0  fc=0 

Comparison  of  coefficients  gives  the  equation 

n 

(n  +  l)6n+1  =  E  (fc  +  l)at+16n_t  (48) 

*=o 

from  which  the  coefficients  can  be  computed  by  recurrence.  Let  1/T(1  +  z)  be  expressed 
by  a  series  in  accordance  with  the  equation 


_ 5 _  =  f  c  zm 

r(l  +  z)  m=o  m 

The  coefficients  are  related  by  the  equation 

cm  =  bm-l  +  bm 

The  first  two  coefficients  are  given  by  the  equations 
c0  —  1 


(49) 


(50) 


(51) 


A  high-precision  value  of  y  was  obtained  from  the  published  literature19  20 
High  precision  values  of  the  coefficients  in  the  convergent  series  are  given  for  the 
gamma  function  in  Table  III.  and  are  given  for  the  digamma  function  in  Table  V 
The  logarithm  of  the  gamma  function  is  expressed  by  the  equation 


log  r( z )  (z  |)logz  z  +  |  log  (2n)  +  E 


B, 


(52) 


m«\  2m(2m  l)zJm  ' 

High- precision  values  of  the  coefficients  in  the  asymptotic  series  are  given  for  the 


i 


I 


gamma  function  in  Table  IV,  and  are  given  for  the  digamma  function  in  Table  VI 

Multiple-Precision  Package 

MP,  or  Multiple-Precision  Package  is  a  co'lection  of  subroutines  which  have  been 
developed  by  Brent24  A  multiple  precision  number  is  represented  by  an  array  of 
floating-point  numbers.  The  first  number  in  the  array  defines  the  sign  of  the  number 
The  next  number  in  the  array  gives  the  exponent  of  the  number.  The  remaining 
numbers  in  the  array  express  the  fraction  of  the  number  Each  successive  number  in 
the  fraction  is  the  coefficient  of  a  progressively  more  negative  power  of  the  base  The 
number  of  coefficients  in  the  composition  of  the  fraction  is  unlimited  The  number  of 
bits  in  the  base  may  have  any  value  not  exceeding  half  the  number  of  bits  in  the 
fraction  of  a  floating  point  number.  Otherwise  multiple-precision  multiplication  would 
be  unduly  cumbersome.  Included  in  the  package  are  routines  for  arithmetic  and  special 
functions. 


APPENDIX  C 


MULTIPLE-PRECISION  NUMBERS 


TABLE  1 

Bernoulli  Numbers 


n  Bn 

0  1 .0000000000000000000000000000000000000000000000000000000000000000  E  000 

1  -5.0000000000000000000000000000000000000000000000000000000000000000  E-01 

2  1 .6666666666666666666666666666666666666666666666666666666666666666  E-0 1 

3  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

4  -3  3333333333333333333333333333333333333333333333333333333333333333  E-02 

5  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

6  2.3809523809523809523809523809523809523809523809523809523809523809  E-02 

7  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

8  -3.3333333333333333333333333333333333333333333333333333333333333333  E-02 

9  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

10  7.5757575757575757575757575757575757575757575757575757575757575757  E-02 

1 1  0  0000000000000000000000000000000000000000000000000000000000000000 E  000 

12  -2.531 135531  135531135531135531 135531135531135531135531  135531135531 1  E-01 

1 3  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

14  1  1 666666666666666666666666666666666666666666666666666666666666666  E  000 

1 5  0 OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOCOOOOOOOOOOOOOOOOOO  E  000 

16  - 7.092 1 568627450980392 1 568627450980392 1 568627450980392 1 568627450980  E  000 

1 7  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

18  5.4971  177944862155388471177944862155388471177944862155388471177944  E  +  01 

1 9  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

20  -5.29 12424242424242424242424242424242424242424242424242424242424242  E-02 

2 1  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

22  6.1921231884057971014492753623188405797101449275362318840579710144  E  +  03 

23  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

24  -8.65802531 135531 135531135531 135531135531 135531 135531 135531 135531 13E+04 

25  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

26  1  4255171666666666666666666666666666666666666666666666666666666666  E  +  06 

27  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

28  -2.729823 10678 1609 1 95402298850574712643678 1609 19540229885057471 2643  E  +  07 

29  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

30  6.01 580873900642368384303868 1 748359 1 677 1 400642368384303868 1 748359 1  E-OB 

3 1  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

32  -  1.5 11 63 15767092 1568627450980392 1568627450980392 1568627450980392 156  E+  10 

33  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

34  4  2961464306 11 6666666666666666666666666666666666666666666666666666  E+  1  1 

35  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

36  -  1.371 16552050883327721590879485616327721590879485616327721590B7948E-  13 

37  0.0000000000000000000000000000000000000000000000000000000000000000  E 000 

38  4.883323 1 8973593 1 6666666666666666666666666666666666666666666666666  E+  1 4 

39  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

40  -  1  92965793419400681486326681448632668144B6326681448632668144863266E+  16 

4  1  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

42  8  4 1 6930475736826 1 500055370985603543743076626799557032 1 15171650055E4  17 

43  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

44  4  033807185405945541307681 1594202898550724637681 1 59420289855072463  E4  19 


TABLE  I 
(Continued) 


n  Bn 

45  O.OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOt’OU  E  000 

46  2.1150748638081991605601453900709219858156028368794326241 134751 773 E  •  21 

47  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

48  -1.2086626522296525934602731193708252531781943546649429002370178840  K  •  23 

49  0.0000000000000000000000000000000000000000000000000000000000000000  1  000 

50  7.50086674607696436685572007575757575757575757575757575757575757':)?  E  -  24 

51  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

52  -5.0387781014810689141378930305220125786163522012578616352201257861  E  •  36 

53  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

54  3.65287764848181233351 10430842971177944862155388471 177944862 1 55388 E  -  28 

55  0.0000000000000000000000000000000000000000000000000000000000000000  E 000 

56  - 2.8498769302450882226269 1 464329 106781 609 1 9540229885057471 2643678 1 6  E  •  30 

57  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

58  2.3865427499683627644645981919219214971751412429378531073446327683  E  •  32 

59  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

60  -2.1399949257225333665810744765191097392674151 161723874574218307692  E  •  34 

61  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

62  2.0500975723478097569921733095672310251666666666666666666666666666  h  •  36 

63  0.0000000000000000000000000000000000000000000000000000000000000000  h  000 

64  -2.0938005911346378409095185290027970184709215686274509803921568627  1  -  38 

65  0.0000000000000000000000000000000000000000000000000000000000000000  1  000 

66  2.27526964884635 1 555964926035276926458 1 46996540588980563023392354 9  E  -  40 

67  0.0000000000000000000000000000000000000000000000000000000000000000 E  000 

68  -2.6257710286239576047303049736158202081449000333333333333333333333  E-  42 

69  O.OOOOOOOOOOOOOOOOUOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO  K  000 

70  3.212508210271803251820479230426496524352194 110616730687153222364  4  L  -  44 

71  0  OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOEOOO 

72  -4  159827816679471091391707449526235893668960301 1346470769224934863  E  •  46 

73  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

74  5.69206954820352800238834562 1912105864448051297 1811 666666b66666666  E  -48 

75  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

76  -8. 2 1636294 197845756922906534686 173330 1 4550892762886003333333333333  E- 50 

77  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

78  1.2502904327166993016732339829702895524177196364448477501 1 15129596  E-  53 

79  O.OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOUOO  E  Oi iO 

80  -2  00 1 5583233248370274925329 1 988 1 3298768724220 1 328259 1 59 1 52074 56  1  37 E  •  55 

81  O.OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOn  E  000 

82  3.3674982915364374233396676903338753016219598947193843672321546104  E-  57 

83  0  OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOU  1  000 

84  -5.94709705031354477186604968440515408405790715651069049904704  31005  E  •  59 

85  0  OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOUOOOOOOOOOOOOOOOOOOOOOOOO  E  000 

86  1  101  191032362797755956413079043769160463051144422314886269994  9716  E  •  63 

87  0  OOOOOOOOOOOOOOJOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO  E  000 

88  2  1355259545253501 18865838501904  106567897329873916346921 1804590304  K  •  64 

89  0  0000000000000000000000000000000000000000000000000000000000000000  E  000 

90  4  332889698664 119241961661 30593792062 18451 3685 118091091 44986557880  E  •  66 


TABLE  I! 


’i 


t 


Riemann  Zeta  Function 

n  t(n) 

1 

2  1 .64  19340668482264364724 15166646025 1892 1 89499012067984377355582293 

3  1  20205690315959428539973816151 14499907649862923404988817922715553 

4  1  08232323371  11381915160036965411679027747509519187269076829762154 

5  1 .036927755143369926331365486457034 1 6805708091950191281 1 9741926779 

6  1.0173430619844491397145179297909205279018174900328535618424086640 

7  1.0083492773819228268397975498497967595998635605652387064172631365 

8  1.0040773561979443393786852385086524652589607906498500203291 102026 

9  1.0020083928260822144178527692324120604856058513948887565485966159 

10  1.0009945751278180853371459589003190170060195315644775172577889946 

11  1.0004941886041 194645587022825264699364686064357582086171 191414361 

12  1.0002460865533080482986379980477396709604160884580034045330499521 

13  1.0001227133475784891467518365263573957142751058955098451367026716 

14  1  0000612481350587048292585451051353337474816961691545494827552022 

15  1.000030588236307U2L  1935517285106450625876279487068581775065699328 

16  1.000015282259408651871732571487636722023237388990471531 1531052035 

17  1 .0000076371976378997622736002935630292130882490902626790953798439 

18  1  0000038 1729326499983985646 164462 19397304546972 1B953331 1431744299 

19  1.0000019082127165539389256569577951013532585711448386302359330467 

20  1.0000009539620338727961 131520386834493459437941874105957500564898 

21  1.0000004769329867878064631 167196043730459664466947849376002074873 

22  1  000000238450502727732990003648 1 8675299493504 1 82 1 77965826984  9603 1 

23  1.0000001192199259653110730677887188823263872549977845198586032257 

24  1.0000000596081890512594796124402079358012275039188373027958642469 

25  1.00000002980350351465228018606370506936601 18447309195433123986813 

26  1.0000000149015548283650412346585066306986288647881678859105474359 

27  1.00000000745071 17898354294919810041706041 194547190318825658299932 

28  1.00000000372533402478845705481920401840242323289305929581 15197693 

29  1  0000000018626597235130490064039099454169480616653304692006657748 

30  1  0000000009313274324196681828717647350212198135679551368161850086 

31  1.0000000004656629065033784072989233251220071062691853369473073729 

32  1.000000000232831 1833676505492001455975940495024829822845303110776 

33  1  0000000001 164155017270051977592973835456309516522471727635932565 

34  1.0000000000582077208790270088924368598910630541731226046172159550 

35  1 .000000000029 1 03850444970996869294252278840464 1 0698 1 9874330322562 

36  1.000000000014551921891041984235929632245318420983808894 1240380691 

37  1  0000000000072759598350574810145208690123380592648509255554661077 

38  1.000000000003637979547378651 1902372363558732735126460283848974699 

39  1.0000000000018189896503070659475848321007300650305893096186640705 

40  1.00000000000090949478402638892825331 18386949087538600009908788285 

41  1.00000000000045474737830421540267991 120294885703390452991 14386280 

42  1 .00000000000022737368458246525 1 522682 1 57797869 1 2 1 382982 1 989 1 58725 

43  1.0000000000001 136868407680227849349104838025906437435902842517998 

44  1.0000000000000568434198762758560927718296752406855305715889938835 

45  1.0000000000000284217097688930185545507370494266207436882653098338 
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E  000 
E  000 
E  000 
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E  000 


TABLE  II 
(Continued) 

n  «T(n) 

46  1. 00000000000001421085482803 160676983430714 17395376786986056339419;  bbv 

47  1 .0000000000000071054273952108527128773544799566000227420435936676  K000 

48  1 .0000000000000035527 1369133711 3673298469534059342992 1 456555030626  K  000 

49  1.0000000000000017763568435791203274733490144002795701555085753269 t  00  1 

50  1.000000000000000888178421093081590309609138639138632560887146464  4  KOOO 

51  1.0000000000000004440892103143813364197770940268121336459603070244  t  OOu 

52  1.00000000000000022204460507980419839993200942046539642366543294  38  E  000 

53  1.0000000000000001 1 10223025141066133720544569921382702483222900442  t'000 

54  1.00000000000000005551 11512484548124372373659050943028167235506166  KOOO 

55  1.0000000000000000277555756213612417258163245385406976898489037436  KOOO 

56  1  0000000000000000 1 387778780972523276283909490650022 1 9077 1 8624686 1  E  000 

57  1.000000000000000006938893904544153697446085326249809274835874 1793  KOOO 

58  1.0000000000000000034694469521659226247442714961093346219504706270  KOOO 

59  1 .000000000000000001 73472347604757657204 B9729699375959074 780 54 4 789  KOOO 

60  1.0000000000000000008673617380119933728342055067342951487907141437  KOOO 

61  1.00000000000000000043368086900206504874970235659062413612547801 16  KOUO 

62  1  000000000000000000216840434499721978501 391 01 68320984576 15740 1040 KOUO 

63  1.00000000000000000010842021724942414063012711 16546138258936474378 KOOO 

64  1.000000000000000000054210 108624566454 109187004043886337 1506342238  KOOO 

65  1.0000000000000000000271050543122346883195462131 194977643 1888728 ! 6 K 000 

66  1  OOOOOOOOOOOOOOOOOOO 13552527 156 101 164581 48523399682692832898 18772 f  000 

67  1 .00000000000000000000677626357804518909799529874 155668620598 12585 KOOO 

68  1. 000000000000000000003388 13 1789020796818085703 1 0045083683403 11 585 i.OOil 

69  1.00000000000000000000169406589450979916540649274712486194030364  17  KOOK 

70  1.0000000000000000000008470329472546998348246992609182167522283864  MA  i 

71  1 . 0000000000000000000004235 164736272B333478622704833 57934 4 088 10971  I  uou 

72  1.00000000000000000000021 1758236813619473184420943981800258694 1761  KOOO 

73  1 .000000000000000000000 1058791184068023385226500 1 53923839847069990  K 000 

74  1  000000000000000000000052939559203398703238 1 39 1 23029 1 850558663756 K  0(H; 

75  1  000000000000000000000026469779601 6985296 11341  1 668420387 1 5592556 1  h  00( > 

76  1 .00000000000000000000001323488980084899080309451 0250944989684  3248  K  000 

77  1 .0000000000000000000000066174449004244040673552453323082200 1 1 7 1 37 !  UOU 

78  1 .0000000000000000000000033087224502 12171 5889469563843 1 4 4048092764  K  OOil 

79  1 .00000000000000000000000 1 654361225 1 06075646229923677 1 8 1 0488297723  K  OOu 

80  1 .000000000000000000000000827 1 806 1 2553034440367 1 1 056 1 674  4  072404009 K.  Oui  ■ 

81  1.00000000000000000000000041359030627651609260093824555081412842. .7  ••  U.'" 

82  1.000000000000000000000000206795153138257670439596791934689504 13  U  '■  .  4i 

83  1 .000000000000000000000000 1 03397576569 1 287099328409559 1 74 48609 1  1  07  ; -"O 

84  1 .000000000000000000000000051698788284564313204101332166355513H9.il.!-  or 

85  1  000000000000000000000000025849394142282142681277617708450222269’  Ki  uu 

86  1.000000000000000000000000012924697071 1410667003811261 183318653093 i  uiu 

87  1  000000000000000000000000006462348535570531803438002161  122 1 670666  f-  O'V 

88  1.000000000000000000000000003231  174267785265386134814  1  18026647  1  :  !•  ..." 

89  1  00000000000000000000000000161558713389263252120601140570522727;:.  ! 

90  1.00000000000000000000000000080779356694631620331587381863408997 


TABLE  III 

Convergent  Series  for  Gamma  Function 


n  cn 

0  1 .0000000000000000000000000000000000000000000000000000000000000000  E  000 

1  0.5772156649015328606065120900824024310421593359399235988057672348  E000 

2  -0.65587807 15202538810770 195 15 14539048 12797663804785843472923624456 E  000 

3  -0.0420026350340952355290039348754298 18711 39450040 1 1 060935220656 1 29  E  000 

4  0.1665386113822914895017007951021052357177815022471743405704689031 E000 

5  -0.0421977345555443367482083012B9187391301652684 1898224863769188732  E  000 

6  -0  0096219715278769735621149216723481989753629422521 130021051388626  E000 

7  0  007218943246663099542395010340446572709904800880238318001094781 1  E000 

8  -0.0011651675918590651 121 139710840183886668093337953840574434075052  E000 

9  -0.0002 15241674 1 1495097281572996305364780647824 1 9233783387503502674  E  000 

10  0.0001280502823881 161861531986263281643233948920996936772149005458  E000 

11  -0.0000201348547807882386556893914210216183822948332979791 152611626E000 

12  -0.000001250493482142670657345359473833092242322655621 1539598 153499  EOOO 

13  0.000001 13302723198 1695882374 12962033074494332400483862 10756542955  EOOO 

14  -0.00000020563384 1 6977607 1 03450 1 54 1 300205728365 1 2579026293379453468  E  000 

1 5  0.0000000061 160951044814 158 178624986828553428672758657 1 97 1 23208673  E  000 

16  0.00000000500200764446922293005566504805999 1 3030446 1 274249448 1 7 1 89  E  000 

17  -0  000000001 181274570487020 144588 1265654365055777387595049325875909 EOOO 

1 8  0  000000000 10434267 11691100510491 5403323 1 2250 1 9 1 400709823 1 258 1 2 1 2 1  E  000 

1 9  0.00000000000778226343990507 12540499373 1 1 360777226068086 1 8 1 3929388  E  000 

20  -0  0000000000036968056 1 8642205708 1 878 1 587808576623657096345 1 360995 1  E  000 

2 1  0  0000000000005 10037028745447597901 548 1 32286323 1 802726886069707632  E  000 

22  -0. 000000000000020583260535665067832224295448552374 1 974609 1 0808 1 0 1 4  E  000 

23  -0.000000000000005348 1225394230179823700173187279399489897154781206  EOOO 

24  0  00000000000000 1226778628238260790 158893846622422428 1654557504563  EOOO 

25  -0  0000000000000001 18125930169745876951376458684229783121 1557291804  EOOO 

26  0.000000000000000001 18669225475160033257977724292867407 10884940796  EOOO 

27  0  0000000000000000014123806553180317815558039475667090370863507503  EOOO 

28  -0  0000000000000000002298745684435370206592478580633699260284505931  EOOO 

29  0 .0000000000000000000 1  7 1440632 1 9273374333839633702672570668 1 265606  S  000 

30  0  000000000000000000000 1 33735 1 730493693 114864  781 395 1 22268022875059  K  000 

3 1  -0.000000000000000000000205423355 176667278932502535 1 355733796682037 1. 000 

32  0  0000000000000000000000273603004860799984483 1 5099043309820 148653 1  E  000 

33  -0  00000000000000000000000 17323564459 1051 663905742845 15647797990697  EOOO 

34  -  0  000000000000000000000000023606 1 90244  9928728734345073542753 1 00792  E  000 

35  0.0000000000000000000000000 1 86498294 1 7 1 72944307 18413161 87866689894  E  000 

36  0  0000000000000000000000000022 1 8095624207 1 972043997 1 69 1 36268603797  E  000 

37  0.000000000000000000000000000 1 29778 1 9749479936688244 14486330594 1 65  E  000 

38  0.00000000000000000000000000000! 1 806974749665284062227454 1 550997 1 5  E 000 

39  -0  00000000000000000000000000000 1 1 2458434927708809029365467426 1 4395  E  000 

40  0  000000000000000000000000000000 1 277085 1751  4086620399020667775 1  1 24  E  000 

4 1  -0  00000000000000000000000000000000739 145116961 514082346 1 2893301085  E  000 

42  0.00000000000000000000000000000000001 13475025755421576095416525946  EOOO 

43  0  00000000000000000000000000000000004639 1 3464 1 05872202994480490795  E  000 

44  -0  0000000000000000000000000000000000053473368 184391 988750774 1 8 1 967  E  000 


1 

■t 


TABLE  III 
(Continued) 


n  cn 

45  0.00000000000000000000000000000000000032079959236 1 335262266 1 237279  E  000 

46  -0.000000000000000000000000000000000000004445829736550756882 1 0 1 5903  E  000 

47  -0.00000000000000000000000000000000000000 1311174518881 9887 1 290 1 0584  E 000 

48  0.000000000000000000000000000000000000000 1 6470333525438 1 38868 1 8259  E  000 

49  -0.0000000000000000000000000000000000000000 1056233 178503581 2 1 860056  E  000 

50  0.0000000000000000000000000000000000000000002678442982643049478354  E  000 

5 1  0.000000000000000000000000000000000000000000024247 1 549485 1 78268967  E  000 

52  -0.0000000000000000000000000000000000000000000037365878345356 1 25540  E 000 

53  0.0000000000000000000000000000000000000000000002628332980940195449  E  000 

54  -0.0000000000000000000000000000000000000000000000092981759953768862  E  000 

55  -0  0000000000000000000000000000000000000000000000002327942418699470  E000 

56  0.0000000000000000000000000000000000000000000000000616962083524438  E  000 

57  -0.0000000000000000000000000000000000000000000000000049282955867709  E  000 

58  0.0000000000000000000000000000000000000000000000000002 1 835 1 3 1 834 1 4  E  000 

59  -0.00000000000000000000000000000000000000000000000000000 1 218722 1 89 1  E  000 

60  -0.0000000000000000000000000000000000000000000000000000007 1171 0884 1  E  000 

6 1  0.0000000000000000000000000000000000000000000000000000000692050405  E  000 

62  -0.0000000000000000000000000000000000000000000000000000000036764384  E  000 

63  0.0000000000000000000000000000000000000000000000000000000000856309  E  000 

64  0.0000000000000000000000000000000000000000000000000000000000049630  E  000 

65  -0.0000000000000000000000000000000000000000000000000000000000007154  E  000 

66  0. 0000000000000000000000000000000000000000000000000000000000000455  E  000 

67  -0.0000000000000000000000000000000000000000000000000000000000000016EOOQ 


ft  -* 


TABLE  IV 


27  1 

28  9 


36  8 


42  8 


Asymptotic  Series  for  Gamma  Function 

_ 5  2m _ 

2m(2m  -  l) 

3333333333333333333333333333333333333333333333333333333333333333  E-02 
7777777777777777777777777777777777777777777777777777777777777777  E-03 
9365079365079365079365079365079365079365079365079365079365079365  E-04 
9523809523809523809523809523809523809523809523809523809523809523  E-04 
4175084175084175084175084175084175084175084175084175084175084175  E-04 
9175269175269175269175269 175269 175269 1752691752691 752691 75269175  E-03 
4102564102564102564 102564 102564 102564 102564 102564 102564 102564102  E-03 
9550653594771241830065359477124183006535947712418300653594771241  E-02 
7964437236883057316493649001588939669435025472177174963552672531 E-01 
3924322169059011164274322169059011 16427432216905901 1164274322169  E  000 
3402864044 1 6839 1994478951 000690 13112491 3733609385783298826777087  E+ 0 1 
56848284626002017306365 13245208897382810426288687 1 58252375643679  E+02 
1931033333333333333333333333333333333333333333333333333333333333 E+03 
6108771253724989357173265219242230736483610046828437633035334184  E  +  04 
9147226885131306710839525077567346755333407168779805042318946657  E+05 
5238221539407416192283364958886780518659076533839342188488298545  E+07 
829007513914 141414141414 14141414 1414 14 14 14 14 14 14 14 14 14 14 14 14 14 14  E+08 
088226603578439 1 0890 151491 65525 1 053747294348798 1 08 1 9660443720594  E+  1 0 
4732028376500225225225225225225225225225225225225225225225225225  E+  1 1 
.2369602 14226927445425 1 7 1 034927 1 32488 1 08097864 195425171 034927 1 324  E+  1 3 
887880647930793350758 1516251 8022902 1 0847053890567382 1 80703629532  E+  1 4 
1320333960919373896975058982136838557465453319851702055948769801  E+ 16 
02 1 7752965257000775652876280535855003940 1 1 032308904649330 1 8 1 2450  E+  1 8 
.3575472173300203610827709191969204484849040543658816499867814010  E+  19 
061578263704883415043151 051329622758194 1867656 1533704390B4724799  E+21 
89999 1 74263992040502937 1 429306942902947342458996 1 77087 1 870760882  E+  23 
.2763374033828834149234951377697825976541633608829901448239746816  E+25 
25284717612041630723024234834762277951933 1 243469 1 745036572622779  E+26 
21 8822595 1 856 1 0297836050 1 8730 1 6379224898404202596887699474675389  E  +  28 
0451 834059958569677431 48238754547286066 1443959671 9620740630 1 6080  E  +  30 
4206704715700945451934778148261000136612021857923497267759562841  E+32 
1 929578153140819467001947643918576846997062713974478680361033302  E+34 
3036588551197005966548392430697586436992926354055490979566255361  E+36 
7633253481649640138944358507809925551907375621890547263681592039  E+38 
651 1557 148484539375 16520 1458105559510397393594549289589093630734  E+40 
.137378358 1366805387 16 172632093575691 8406891649738792623679450044 E+42 
0536966953357141803754804927641810189648373375011415525114155251  E+45 
4418180599962206261805377801511812809570332063664211 111111111111  E+47 
081735652208956546242480824126356231 1317343264149979189335880 1 13  E+49 
167022663488666 18274 1349556774256 1 34291 80698304207530303915446 16  E+51 
.0700064612111373431792648153174876567629628044555621307319401061  E  +  53 
52997282030055 1 8816208400522 1 62278887807044700382824 1 40089545446  E  +  55 
.50641 72809340598576695 1  1 7360379879076 1 0 1 93084 02505456398084 1 2745  E+  58 
789349470383 163687128838 16863 1278 171 2347568886054688 103193038537  E+60 
40935043528604 1 500576356 18718841 525823362902769 1 524237826037 1 885  E+  62 


TABLE  V 

Convergent  Series  for  Digamma  Function 
n  ncn 

0  0.0000000000000000000000000000000000000000000000000000000000000000  E  000 

1  0.57721566490 153286060651209008240243 1042 1593359399235988057672348  E  000 

2  -1.31 1756143040507762 154039030290780962559532760957168694584724891 3 E 000 

3  -0.12600790510228570658701 18046262894561341 8350120331 82805661 974389  E  000 

4  0.666154445529165958006803180408420942871 1260089886973622818756127 E000 

5  -0.2109886727777216837410415064459369565082634209491 124318845943663  E000 

6  -0.0577318291672618413726895300340891938521776535126780126308331757  E  000 

7  0.05053260272664 1696796765072383126008969333606161 6682260076634682 E  000 

8  -0.00932 134073487252089691 176867214710933447467036307245954 72600422  E  000 

9  -0.00193717506703455875534 1569667482830258304 177310405048753 1524074  E  000 

10  0.001280502823881 1618615319862632816432339489209969367721490054583 E000 

11  -0.0002214834025886706252125833056312400022052431662777702678727893 E  000 

12  -0.00001500592 1785712047888 1443 13685997 106907B7 1 86745384751 7784 1 990  E  000 

13  0.00001472935401576204647086368506429968426321206290207398350584 15  E000 

14  -0.00000287887378376864994483021 5782028801 97 1 1 1 76106368 1 073 1 2348556  E 000 

15  0  0000000917414265672212372679374802428301430091379857956848130098  E  000 

16  0.00000008003212231150756688089064076895986084871380387991 1 7075032 E  000 

17  -0.00000002008166769827934245799815161242059482155891 15838539890463  E000 

18  0.00000000 18781680810439809 1888477259816205034452 127768 1626461 81 79 E 000 

19  0.000000000 14786300535819635382694880891 58547672952936374  464658375 E  000 

20  -0.000000000073936 1 1 2372844 1 141637563 1756171 532473 1 4 1 92690272 1 99027  E  000 

21  0  0000000000107107776036543995559325107780127667857264607463860274  E000 

22  -0.0000000000004528317317846314923089344998681 52232344 1400377782323  E  000 

23  -0.000000000000 1 230068 1 84067294 1 35945 1 03983307426 1 88267634559967756  E  000 

24  0.0000000000000294426870777 182589638 1 34523 1 8938 1 382759709380 1 095 1 7  E  000 

25  -0.00000000000000295314825424364692378441 14671057445780288932295121  E000 

26  0  0000000000000000308539986235416086470742083 1 6 1 4552584830081607 12  E  000 

27  0  000000000000000038 1 34277693586858 1 0200670658430 1 14400 1 33 14702590  E  000 

28  -0.0000000000000000064364879 1641 90365784589400257743579287966 1 66079  E  00 0 

29  0.0000000000000000004971778333589278556813493773775045493756702581  E000 

30  0  0000000000000000000040 1 2055 191481 079344594344 1 85366804068625 1 784  E  000 

31  -0.000000000000000000006368 1240 104766856469075785892027747697 1 43 1 75  E  000 

32  0.0000000000000000000008755296 1 5554559950346083 1 693859 1 42447568997  E  000 

33  -0.000000000000000000000057 1 677627 1 5047049088895 1 3890 163773336930 1 7  E  000 

34  -0.0000000000000000000000008026 104683297576776967732500453605426949  E  000 

35  0.00000000000000000000000065274402960 1 05305075 14446066575334 1 463 1 0  E  000 

36  -0  0000000000000000000000000798514424714590993583898088905669736703EOOO 

37  0.0000000000000000000000000048017933073075765746503335994  231984 127  E  000 

38  0.0000000000000000000000000000448665040487280794  364643257893789 1 77  E  000 

39  -0  000000000000000000000000000043858789621 80643552 1 452532296 1 96 1 409 E 000 

40  0  00000000000000000000000000000510834070056346481596082671 1 004  4985  E  000 

4 1  -0.00000000000000000000000000000030304949795422077376 1 9 1 28625344506  E  000 

42  0  000000000000000000000000000000000476595 1 08 1 727706 1 960074  9408977 1  E  000 

43  0.0000000000000000000000000000000019948278956552504728762661 104194  E  000 

44  -0.00000000000000000000000000000000023528282001  1 3247505034064006552  E  000 


TABLE  V 
(Continued) 

n  nc„ 

45  0.0000000000000000000000000000000000 1 443598 1 6562600868028755677558  E  000 

46  -0.000000000000000000000000000000000000204508 1 6788 1 334816576673 1 56 1  E  000 

47  -0.0000000000000000000000000000000000000616252023874534695063497492  E  000 

48  0.0000000000000000000000000000000000000079057600922 103066567276447  E  000 

49  -0.0000000000000000000000000000000000000005 17554257466754797 1 142749  E  000 

50  0.0000000000000000000000000000000000000000 1 33922 1491321 524739 1 7748  E  000 

51  0.0000000000000000000000000000000000000000012366049023744091 717332 E  000 

52  -0.000000000000000000000000000000000000000000 1 9430256739585 1 8528097  E  000 

53  0.0000000000000000000000000000000000000000000 1 3930 1 64798983035880 1  E  000 

54  -0.0000000000000000000000000000000000000000000005021015037503518601  E000 

55  -0.0000000000000000000000000000000000000000000000128036833028470882  E  000 

56  0.0000000000000000000000000000000000000000000000034549876677368569  E  000 

57  -0.0000000000000000000000000000000000000000000000002809128484459464  E  000 

58  0.0000000000000000000000000000000000000000000000000 1 2664376463804 1  E  000 

59  -0.000000000000000000000000000000000000000000000000000071904609 1597  E  000 

60  -0.0000000000000000000000000000000000000000000000000000427026530499  E  000 

61  0.0000000000000000000000000000000000000000000000000000042215074731  E000 

62  -0.0000000000000000000000000000000000000000000000000000002279391850  E  000 

63  0.00000000000000000000000000000000000000000000000000000000539475 1 7  E  000 

64  0.0000000000000000000000000000000000000000000000000000000003 1 76349  E  000 

65  -0.0000000000000000000000000000000000000000000000000000000000465029  E  000 

66  0.000000000000000000000000000000000000000000000000000000000003004 1  E  000 

67  -0.000000000000000000000000000000000000000000000000000000000000 1 084  E  000 

68  -0.0000000000000000000000000000000000000000000000000000000000000002  E  000 

69  0.0000000000000000000000000000000000000000000000000000000000000003  E  000 
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Asymptotic  Series  for  Digamma  Function 
2m 

8  3333333333333333333333333333333333333333333333333333333333333333 E-02 

8.3333333333333333333333333333333333333333333333333333333333333333  E-03 
3. 9682539682539682539682539682539682539682539682539682539682539682  E-03 

4. 1666666666666666666666666666666666666666666666666666666666666666  E-03 
7.5757575757575757575757575757575757575757575757575757575757575757  E-03 
2. 1092796092796092796092796092796092796092796092796092796092796092  E-02 

8.3333333333333333333333333333333333333333333333333333333333333333  E-02 
4.4325980392 1 568627450980392 1 568627450980392 1 568627450980392 1 56862  E-  0 1 
3.0539543302701 197438039543302701 1 97438039543302701 1 97438039543302  E  000 
2.6456212121212121212121212121212121212121212121212121212121212121  E+01 
2  8146014492753623188405797101449275362318840579710144927536231884  E  +  02 
3.6075105463980463980463980463980463980463980463980463980463980463  E+03 

5.4827583333333333333333333333333333333333333333333333333333333333  E+04 
9.749368238505747 1 2643678 1 609 1 954022988505747 1 2643678 1 609 1 95402298  E+  05 
2  0052695796688078946 1 4346227249453055904668B078946 1 43462272494530  E+ 07 
4.7238486772 1 62990 1 9607843 1 372549019607843 1 3725490 1 96078431 3725490  E  +  08 

1 .26357247959 1 6666666666666666666666666666666666666666666666666666  E+  1 0 
3.8087931 12524536881 155302207933786881 1553022079337868811553022079  E+  1 1 

1 .2850850499305083333333333333333333333333333333333333333333333333  E+  13 
4.824 1448354850170371 581670362158167036215816703621581670362158167  E+ 14 
2.00403 1 06565 1 625273810842 1 663238938986447292095 1 32626694088488 1 08  E+  1 6 
9. 167743603 19533077569927536231884057971014492753623188405797 10144  E+ 18 
4.5979888343656503490437943262411347517730496453900709219858156028  E+ 19 
2.5160471921451 095697089023320225526 1 078790490555 1 9643754937672584  E+  2 1 

1  5001733492153928733711440151515151515151515151515151515151515151  E+23 

9  6899578874635940656497942894654088050314465408805031446540880503  E+24 
6.76458823792928209909452423017984776756706581267984776756706581 26  E+26 
5. 08906594686622896897663329 159 1 1925287356321 8390804597701 14942528  E  +  28 
4.1147288792557978697665486067619336158192090395480225988700564971  E  +  30 
3.5666582095375556 10968457460865 1 82898779025 1 936206457623697 1 79487  E  +  32 
3.30660898765775767256802 146704392 10083333333333333333333333333333  E+34 
3.271 56342364787 1 62642 1122701 56687034 1 3608 1 4950960392 1 568627450980  E+  36 
3.447378255827805387825645507995343 1 1 84045402 1 30136069 1 36718065984  E-f  38 
3.86 1 4279832705258893092720200232650 1 1 9777941666666666666666666666  E4  40 
4 .5892974432454332 1 68863989006092836062 174201 5802390098 1 6474605206  E+  42 
5  7775386342770431824884825687864387412068893071314542762812409531  E+44 
7  6919858759507135167410075971785214384433125637583333333333333333  E+46 
1 .08 1 363544997 1 65469635403335 1 1 33859607 1 77749047748 1 58333333333333  E+  49 
1 .60293645220089654060671023457729429797 1 43543 1 3395483975788627687  E  +  5 1 
2.501 947904 156046284365666 1 498516623460905275166032394894009320246  E+53 
4  1 0670523358 1 02 1 247975204500407 165001 97799987 1 60900532589287 1 4859  E+  55 
7 07987744084945806 1 745297243339469 1 476879847 1 0 1 3 1 7744036274322720  E+  57 
1 .2804546887939508790 1 908497563228972 1 468664 121421 2963793837 1 50833  E+60 

2  4267340392333524078020892067092 120089742384930867578649777943527  E  +  62 
4  8 1 432 1 8874045769355 1 295700659768957982792983464545657 1 665 1 730978  E  +  64 


1 


1 


APPENDIX  D 


DISTRIBUTION 


i 


DISTRIBUTION 

Defense  Documentation  Center 
Cameron  Station 
Alexandria,  VA  22314 

Defense  Printing  Service 
Washington  Navy  Yard 
Washington,  DC  20374 

Library  of  Congress 

Washington,  DC  20540 

Attn:  Gift  and  Exchange  Division 

C 

D 

E 

F 

G 

K 

K05H 

K05S 

K05D 

K10 

K20 

K30 

K40 

K50 

K60 

K70 

K80 

R 

U 

V 

X21 

X220 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  OF  THIS  PACE  (When  Dale  Entered)  _ 

REPORT  DOCUMENTATION  PAGE  befo^compTenTform 

r  REPORT  NUMBER  {2.  GOVT  ACCESSION  NO.  3  BECIPlENT'S  CATALOG  NUMBER 


|  I.  REPORT  NUMBER 

NSWC/DL  TR-3788 


4.  TITLE  (and  Subtitle) 


5.  TYPE  OF  REPORT  4  PERIOD  COVERED 


COMPUTATION  OF  SPECIAL  FUNCTIONS 


6.  PERFORMING  ORG.  REPORT  NUMBER 


7.  authors; 


8.  CONTRACT  OR  GRANT  NUM8E Rf tj 


A.  V.  Hershey 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Naval  Surface  Weapons  Center  (K05) 
Dahlgren,  Virginia  22448 


II.  CONTROLLING  OFFICE  NAME  ANO  ADDRESS  '2.  REPORT  DATE 

November  1978 

13.  NUMBER  OF  PAGES 

76 _ 


4.  MONITORING  AGENCY  name  4  AODRESSri/  different  from  Confrollfn*  Office)  IS.  SECURITY  CLASS,  fof  Ihlm  roporf) 

UNCLASSIFIED 

'  IS*.  DECLASSIFICATION  DOWNGRADING- 
SCHEDULE 


16.  Distribution  statement  (of  this  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (of  r»i»  abstract  entered  In  Block  20,  II  dlllarant  from  Report) 


19.  KEY  WORDS  7 Continue  on  ravaraa  aide  It  neceeeary  and  Identify  by  block  number) 


potentials 

disk 

plate 

pgendre  functions 


Bessel  functions 
analysis 
programming 
computation 


20  ABSTRACT  fConflnu*  on  ravaraa  aide  If  nocwoawy  and  Identify  by  bfocA  number) 

•^Documentation  is  given  for  some  subroutines  which  compute  potentials  and 
other  functions.  A  set  of  subroutines  uses  rational  approximations  to 
compute  Bessel  functions  of  integral  order.  One  subroutine  uses  the  Debye 
approximation  for  the  efficient  computation  of  Bessel  functions  of  complex 
argument  and  complex  order.  Kmpirical  formulae  have  been  developed  to 
express  the  limiting  boundaries  of  the  modes  of  computation. 


DO  i"RM73  1473  EDITION  OF  1  NOV  65  IS  OBSOLETE 

S'N  0102-LF-014-6401 


UNCLASSIFIED 

s E CURlTY  CLASSIFICATION  OF  THIS  PAGE  (When  Date  tnlered) 


